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Summary 


Each  item  of  work  described  in  this  report  was  paid  for,  in  whole  or  in  part, 
by  the  present  contract.  There  are  four  items.  The  first  item  is  about  nonlinear  mi¬ 
crowave  amplification.  The  second  item  is  about  speeds  of  frequency  components  of 

electromagnetic  pulses,  and  a  related  inversion  method.  The  third  item  is  about  an 

« 

inversion  method  for  depth-dependent  dispersive  media,  which  uses  electromagnetic 
pulses.  The  fourth  item  is  about  general  rules  of  pulse  propagation  in  dispersive 
media,  with  applications  to  inversion. 

Some  of  the  work  was  done  jointly  with  others  or  with  some  help  from  others. 
I  will  now  list  these  people  and  identify  those  who  work  in  San  Antonio.  For  the  first 
item  I  would  like  to  thank  Gary  Thomas  who  told  me  about  an  error  in  an  equation 
and  its  correction.  A  project  with  the  Naval  Weapons  Support  Center  required 
rederivations  done  in  this  first  work  item.  The  second  item  was  done  jointly  with 
Mike  Hobart.  Some  computer  programs  came  from  Robert  Krueger  and  Ronald 
Winther,  and  I  am  grateful  for  conversations  about  error  bars  with  Phelps  Crump, 
who  is  a  defense  contractor  working  in  San  Antonio.  The  fourth  item  was  done 
jointly  with  Peter  Petropoulos  while  he  was  employed  by  Optech,  and  Fred  German 
allowed  us  to  use  his  finite-difference  program. 

This  summary  concludes  with  a  brief,  item-by-item  description  of  the  re¬ 
mainder  of  this  report. 
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Comment  on  a  Paper’s  Nonlinear-Amplification  Model 
The  paper  commented  on  is  significant  because  its  claims  of  laboratory  success 
have  motivated  further  use  of  nonlinear  Schrodinger  equations  to  model  amplifiers. 
Some  of  the  paper’s  inaccurate  equations  are  corrected  here,  and  others  identified. 
Questions  are  also  reused  about  the  claims  of  success  in  the  laboratory. 

Speeds  and  Inversion 

Our  numerical  computation  of  precursors  shows  that  damping  increases  the  speed 
of  frequency  components  to  a  degree  that  is  consistent  with  energy  velocity.  The 
speed-up  here  is  slight,  and  a  medium  that  has  been  called  “highly  absorptive” 
exhibits  no  speed-up.  These  small  effects  lead  us  to  question  the  extent  to  which 
energy  velocity  would  be  more  useful  than  group  velocity  in  the  laboratory.  We  also 
introduce  an  especially  simple  method  for  measuring  permittivity.  The  method  is 
useful  when  damping  affects  wave  speed  no  more  than  slightly,  as  happens  in  all 
numerical  simulations  cited  or  reported  here. 

Inversion  Method  for  Depth-dependent  Dispersive  Media 
This  item  describes  an  inverse  scattering  algorithm  for  electromagnetically  disper¬ 
sive  objects  that  are  flat,  and  whose  properties  vary  only  with  depth.  Tissue  and 
soil  are  discretely  layered  examples  of  such  media,  but  this  paper’s  algorithm  edso 
applies  to  dispersive  media  that  are  continuously  layered.  The  inversion  algorithm 
uses  time-dependent  reflection  data  as  a  function  of  the  angle  of  incidence.  These 
data  are  represented  by  a  function  of  two  variables  A  count  of  the  variables 


2 


suggests  that  the  data  R^{t)  should  be  sufficient  to  determine  the  two- variable 
function  g{x^t)  that  represents  depth-dependent  dispersion.  A  time-domain  layer 
stripping  method  that  uses  wave  splitting  and  Krueger-Ochs  Green  functions  is 
presented. 

Propagation  and  Related  Topics,  Including  Inversion 

We  develop  some  new  methods  for  describing  pulse  propagation  for  general  disper¬ 
sive  media,  using  a  Debye  model  for  water  as  an  example.  Short-pulse,  long-pulse, 
short-time,  and  long-time  approximations  are  presented.  We  explain  a  factor-of- 

i 

nine  effect  in  the  speed  of  waves  in  water,  which  seems  to  have  been  previously 
unnoticed.  We  ^dso  study  the  following  problem:  Knowing  only  the  peak  amplitude 
and  power  density  of  an  incident  pulse,  what  can  be  said  about  the  peak  amplitude 
of  the  propagated  pulse?  We  provide  sharp  upper  bounds  for  the  propagated  am¬ 
plitude  and  reduce  the  computation  of  those  bounds  to  a  Ccdculator  exercise.  These 
bounds  may  be  useful  in  controlling  the  electromagnetic  interference  or  damage 
produced  in  dispersive  media. 
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Comment  on  a  Paper’s  Nonlinear- Amplification  Model 


f 


Continued  citation  of  Gary  E.  Thomas’s  paper  [1]  makes  comment  worth¬ 
while.  In  this  comment,  the  numbers  of  equations,  figures,  tables,  and  sections  are 

of  those  in  [1]. 

Equations 

This  section  identifies  twelve  Inaccurate  equations  in  [1],  which  are  corrected 
here  or  proven  inaccurate.  The  inaccurate  equations  are  (10a),  (10b),  (11)— (13), 
(15a),  (15b),  (21),  the  first  two  unnumbered  equations  that  follow  (15b),  and  two 
of  the  equations  between  (33)  and  (34).  The  signs  of  the  paper  s  functions  J^j,  —Bq, 
and  -Co  aifect  the  accuracy  of  (34)  and  the  equation  for  r  that  is  between  (33) 


and  (34).  Those  two  equations  are  replaced  here  by  ones  that  are  valid  regardless 
of  the  signs  of  Ji,  —Bq,  and  —Co- 

Each  equation  in  the  following  list  corrects  a  corresponding  equation  in  [1], 
or  is  used  immediately  below  to  correct  a  larger  concept. 


E  =  j4sech[A(x  —  u^t)lV2\  exp[iUe(x  —  Uct)l‘l\  (10a) 

E  =  Btanh[B(x  +  Uct)/V2]  exp[iuc{x  +  Uct)l2]  (10^) 

dlf  +  6(sech20/  -  (1  -  2^|M")/  ~  +275M'  (15a) 

dig  +  2{sech^09  -  (1  -  ‘2:kl/A^)g  ^  -2-ff/A^  (156) 

^  =  A{x  —  Uet)  /  V2  {15c) 

=  2kl/A^  and  T  =  27/ (15d) 

tte  =  Uc  i  y/ Uc^  +  2A^  (21) 

iadtE  +  bdlE+cdlE  +  J\E\'^E  =  0  (33a) 

a  =  — 2<jjDq  b  —  — 2Bq1^  c  —  Cq  J  —  J\  (336) 

T  =  t\J/a\  ^  =  xy/\J/b\  ^  =  zy/\jjc\  (33c) 

sign(a/J)  =  1  (33d) 

idrE  +  [sign{b/J)]dlE+[s\sn{c/J)]dlE  +  \E\^E  =  0  (34) 


In  (10a):  A,  Ue,  and  Uc  are  real- valued. constants  that  satisfy  Ue(2uc  —  u^)  = 
2A^  and  u^{ue  —  2uc)  >  0.  In  (10b):  B,  u^,  and  Uc  are  new  real-valued  constants 
that  satisfy  Uf[2uc  —  Ue)  =  2B^  and  u^[2uc  —  u^)  >  0. 
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Table  I  and  Fig.  3  are  from  [1,  Ref.  8] — which  is  the  eighth  reference  in  [l] — 
provided  (15d)  replaces  similar  equations  in  [l],  and  provided  (1,  Ref.  8]  uses,  in  its 
notation,  iV^  =  1.  Reference  [1]  does  not  say  whether  [1,  Ref.  8]  uses  =  1. 

Equations  {33)-(34)  correct  a  corresponding  set  of  equations  in  [1]  only  if  a, 
6,  c,  and  J  are  real  valued  and  independent  of  the  real  variables  t,  i,  and  z.  Equation 
(33d)  is  not  in  [1];  the  equation  comes  from  [2].  There  are  no  more  corrections. 

Equation  (11)  is  false  because  it  and  definition  (9)  have  identical  right  sides, 
though  the  left  sides  have  different  units.  Equation  (11)  is  not  from  (1,  Ref.  7],  nor 
does  (11),  (12),  or  (13)  follow  from  direct  application  of  [1,  Ref.  7].  The  right  sides 
of  (9),  (12),  and  (13)  have  identical  units,  though  the  units  of  the  left  sides  of  (12) 
and  (13)  differ  from  the  units  of  the  left  side  of  definition  (9).  That  is,  (11)— (13) 
have  errors.  Those  equations  are  a  foundation  for  the  paper’s  predictions  about 
power. 

Some  work  reported  after  publication  of  [1]  relates  to  foundational  equations 
mentioned  in  the  previous  paragraph.  In  particular,  much  of  (8)-(13)  concerns 
the  energy-like  term  (9),  which  is  said  to  be  conserved  for  a  nonlinear  Schrodinger 
equation  with  one  spacial  variable,  and  to  be  unconserved  for  spacial  dimension 
two.  Later  work  [3]  shows  that  idtE  +  d^E  +  dlE  +  \E\^E  =  0  implies  that  the 
different  energy-like  term  JJdxdz[\dxE\^  +  \dzE\^  —  |jE|^/2)  is  conserved  and  the 
term  JJdxdz\E\^^  which  is  proportional  to  an  integral  of  the  right  side  of  (13),  also 
is  coriserved. 

This  section  has  identified  twelve  inaccurate  equations  and  corrected  most  of 
them.  The  focus  has  been  on  equations  themselves,  not  their  relevance  in  modeling 
amplifiers.  The  strongest  argument  in  [l]  for  the  relevance  of  its  equations  is  that 


they  are  said  to  make  predictions  that  are  validated  experimentally.  However, 
the  degree  to  which  this  comment’s  corrections  would  affect  those  predictions  is 
unknown  because  [l]  omits  numerical  values  of  some  important  variables.  (The 
variables  are  used  in  (36)-(38)  and  (44)  to  define  apparently  significant  signs  of 
coefficients  in  a  nonlinear  Schrodinger  equation  (33).)  That  is,  predictions  in  [1]  are 

questionable  because  of  errors,  but  the  magnitude  of  the  errors  is  unknown  because 

« 

the  predictions  are  not  reproducible. 

Experimental  Validation 

The  model  in  [1]  was  derived  nonrigorously,  so  its  scientific  relevance  is  based 
on  clziims  of  experimental  validation.  The  paper  says  that  predictions  and  experi¬ 
ments  agree  closely,  so  if  the  numerical  effect  of  errors  in  [l]  is  large,  then  there  is  no 
need  to  examine  [1]  further.  However,  if  the  errors  affect  predictions  only  slightly, 
then  the  scientific  relevance  of  [1]  relies  on  claims  that  (flawed)  predictions  agree 
with  experiments.  This  section  examines  ciU  such  cl2Lims  of  experimental  validation. 

Figure  10  and  Tables  II  and  III  contain  all  the  experimental  data  in  [1], 
and  Table  II  is  illustrated  in  Figs.  6  and  7.  Every  claim  of  experimental  validation 
naturally  refers  to  those  data.  The  claims  are  in  the  second  and  third  paragraphs 
of  Section  V,  and  in  the  fourth  and  fifth  paragraphs  of  Section  VI.  We  will  examine 
claims  in  that  order. 

The  claims  of  experimental  validation  that  concern  Tables  II  and  III  are 
affected  by  “transition  shifts.”  The  shifts  change  numerical  predictions.  Shifts  are 
explained  qualitatively  in  [1],  but  their  numerical  values  are  unexplained.  One  shift 
in  Fig.  4  changes  a  prediction  by  an  unexplained  32%. 
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Let  us  examine  at  length  a  single  prediction.  Section  V  predicts  an  increase 
in  the  slope  of  a  graph  at  the  voltage  0.63Vh,  which  is  shifted  by  17%  to  0.74Vh 
in  Fig.  4.  The  text  says  that  experimental  data  have  an  increasing  slope  at  the 
nearby  voltage  0.79Vh,  and  Fig.  7  is  cited.  Linear  interpolation  is  used  in  every 
interval  of  Fig.  7,  except  for  one  that  has  an  unexplained,  nonlinear  interpolation. 
The  initial  slope  of  the  nonlinear  interpolation  is  chosen  to  agree  with  the  predicted 
increase;  yet  the  slope  decreases  to  the  right  of  0.79Vhj  and  throughout  the  interval. 
Fig.  6,  which  has  the  same  data  as  Fig.  7,  illustrates  the  decreasing  slope  more 
plainly.  Thus,  two  graphs  contradict  a  predicted  increase  in  slope,  while  the  text 
claims  close  agreement.  In  another  matter  concerning  the  same  prediction,  a  table 
does  not  distinguish  experiment  from  theory.  In  particular,  the  voltage  0.79Vh 
comes  from  analysis  of  experimental  data,  but  Table  II  portrays  that  voltage  as  a 
prediction.  A  larger  set  of  related  data  cdso  does  not  clearly  establish  that  theory 
and  experiment  agree.  For  instance,  a  7.2kV  “predicted  value”  is  placed  in  the 
same  row  of  Table  II  as  a  7.2kV  applied  voltage,  a  —  13.7dBm  power,  a  240.0  degree 
phase  shift,  and  a  0.02A  current;  but  the  relation  of  the  7.2kV  prediction  to  any  of 
the  four  measurements  in  its  row  is  unexplained.  Is  the  unexplained  placement  of 
a  7.2kV  prediction  meant  to  support  the  text’s  claim  that  theory  and  experiment 
agree?  Does  that  placement  rely  on  the  extraordinary  interpolation  technique  that 
earlier  produced  two  contradictions?  Those  questions  can  be  repeated  for  the  other 
predictions  in  Table  II,  which  are  also  affected  by  shifts.  Questions  have  now  been 
raised  about  the  claims  of  experimental  validation  that  refer  to  Table  II.  Those 
claims  are  in  the  second  paragraph  of  Section  V. 


Only  one  claim  of  experimental  validation  [1,  page  32]  mentions  Table  III. 
The  claim  refers  to  a  voltage  1.13Vhj  which  is  shifted  twice:  first  by  8%  to  1.22Vh; 
second,  by  an  unmentioned  quantity  whose  sign  might  be  implied  by  the  text, 
which  says,  “For  a  forward-wave  CFA,  negative  phase  velocity  effects  increase  with 
increasing  frequency.”  This  doubly- shifted,  indeterminate  voltage  is  said  to  explziin 

a  feature  in  Table  III.  The  resulting  claim  of  experimented  validation  is  unconvincing 

« 

because  the  numerical  value  of  the  first  shift  is  unexplained  and  the  numerical  value 
of  the  second  shift  is  unmentioned. 

Questions  have  now  been  raised  about  the  claims  of  experimental  Vedidation 
that  concern  Tables  II  and  III.  All  validation  claims  wiU  have  been  examined  once 
we  consider  those  that  refer  to  Fig.  10.  The  claims  about  Fig.  10  are  in  the  fourth 
and  fifth  paragraphs  of  Section  V.  The  subject  is  phzise  shifts. 

Phase  predictions  in  (23)  are  affected  by  an  error  in  (21),  which  is  corrected  in 
the  previous  section  of  this  comment.  It  is  not  clear  whether  corrected  predictions 
would  agree  with  the  data  presented.  The  paper  says,  “phase  shifts  are  scaled 
according  to  power  output,”  but  the  paper  does  not  further  define  a  scaling  formula 
with  which  the  variable  A  in  (21)  could  be  replaced  by  the  power  graphed  in  Fig.  6. 
A  scaling  formula  would  cdlow  the  phase  in  (23)  to  be  graphed,  but  the  paper  offers 
only  qualitative  analysis  and  sketches.  Scaling,  which  is  not  fully  explained,  is  a 
crucial  part  of  claims  that  qualitative  anadysis  and  data  agree  closely.  Here  and 
in  the  case  of  transition  shifts,  experimental  validation  relies  on  the  accuracy  of 
incompletely-reported  work. 
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Conclusion 


Reference  [1]  apparently  is  the  first  paper  to  use  a  nonlinear  Schrodinger 
equation  to  describe  microwave  amplifiers.  The  model  in  [1]  is  is  not  rigorously  de¬ 
rived;  instead,  the  model’s  relation  to  amplifiers  is  based  on  claims  that  predictions 
and  experiments  agree  closely.  Those  claims  have  motivated  other  papers  that  use 

similcir  nonlinear  models. 

« 

This  comment  has  examined  some  equations  that  underly  predictions,  and 
it  hzis  reviewed  all  claims  of  experimental  validation.  Questions  have  been  raised 
about  both  the  equations  and  the  claims  of  vailidation. 

Although  no  attempt  has  been  made  to  re-derive  every  equation  in  [1],  most 
of  the  twelve  known  inaccurate  equations  in  (l]  have  been  corrected  here.  The 
effect  of  corrections  on  numerical  predictions  in  [1]  is  unknown  because  the  values 
of  some  critical  parameters  are  omitted  from  [1].  Another  fundamental  concern  is 
the  paper’s  development  of  an  energy-like  term  (9).  That  development  is  flawed 
and  uncorrected,  yet  it  is  a  foundation  for  predictions  about  an  amplifier’s  power. 

The  effect  of  errors  on  numerical  predictions  in  [1]  has  unknown  magnitude, 
so  it  is  possible  that  the  effect  is  small.  In  that  case,  claims  of  experimented  valida¬ 
tion  are  crucial  because  they  are  the  only  connection  between  [l]  and  the  amplifiers 
it  models.  Most  of  those  claims  involve  “transition  shifts,”  which  are  large  changes 
in  predictions.  The  numerical  values  of  shifts  are  unexplained,  but  not  every  pre¬ 
diction  is  shifted:  zJl  unShifted  predictions  are  “seded”  in  a  way  that  is  also  not 
fully  explained.  Shifts  and  scalings  are  both  documented  incompletely,  and  they 
both  bring  predictions  into  closer  agreement  with  experiments.  The  paper’s  most 
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prominent  claim— the  subject  of  Fig.  7— relies  on  both  an  unexplained  17%  shift 
and  am  extraordinarily  questionable  interpolation. 

In  summary,  several  papers  have  been  motivated  by  the  development  in  [1] 
of  a  nonlinear-Schrodinger-equation  (NLS)  model.  If  the  numerical  effect  of  twelve 
inaccurate  equations  in  [1]  is  not  too  large,  then  the  argument  for  an  NLS  model 
relies  on  shifts,  scalings,  and  questionable  interpolations.  The  shifts  and  scalings 
'  will  not  be  reproducible  unless  [1]  is  further  documented. 
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Speeds  and  Inversion 


1.  INTRODUCTION 

There  are  several  velocities  of  light, among  which  group  and  energy  velocity 
predict  most  successfully  the  speeds  of  a  pulse’s  frequency  components.  This 
paper  compares  those  group-  and  energy-velocity  predictions,  and  says  that  two 
predictions  differ  significantly  only  if  there  is  numerical  evidence  that  the  difference 
could  be  noticed  in  a  laboratory.  A  simple  method  for  measuring  permittivity  is 
•  also  examined  from  in  empiricsd  point  of  view,  and  is  found  to  be  useful  in  some 
cases. 

Group  velocity  =  dw/dk  (with  k  =  u)^/e/c)  is  often  used  for  undamped 
media,  but  is  widely  considered  useless  for  damped  media. First,  we  assert 
to  the  contrary  that  waves  computed^  for  a  ^highly  absorptive^  medium  and 
reported  in  this  journal  are  described  with  equal  accuracy  by  group  and  energy 
velocity.  That  is,  damping  negligibly  affects  the  speed  of  frequency  components  in 
Ref.  5.  Second,  we  report  simulations  in  which  damping  increases  speed  noticeably, 
though  slightly.  Our  experience  of  finding  several  media  in  which  damping  negligibly 
affects  speed  (before  finding  the  exceptional  media  reported  here)  Is  more  anecdotal 
evidence  that  the  effect  of  damping  on  wave  speed  is  slight.  Third,  we  observe 
that  precursors  in  waves  computed  for  a  Debye  medium  and  reported  in  this  journal^ 
are  also  consistent  with  a  hypothesis  that  damping  affects  wave  speed  slightly. 

We  view  those  three  observations  as  the  first  numerical  tests  of  the  empirical 
usefulness  of  energy  velocity.  The  tests  are  reviewed  after  equation  (3),  and  in 
Sections  4  and  5.  We  conclude  that  there  is  no  numerical  evidence  that  energy 
velocity  could  be  much  better  than  group  velocity  at  predicting  speeds  of  frequency 
components  that  are  observable  in  a  laboratory.  However,  the  available  numerical 
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tests  are  not  exhaustive,  and  the  tests  reported  here  tend  to  improve  the  stature 
of  energy  velocity. 

Finally,  we  use  wave  speed  in  an  elementary  technique  for  inferring  permittivity 
from  transmitted  waves.  The  technique  is  truly  elementary:  It  is  not  mathematically 
advanced;  indeed,  we  prove  limitations.  The  elementary  technique  is  worthwhile 
because  it  is  simple,  and  potentially  useful  in  some  cases. 

This  paper  is  mainly  about  Lorentz  media,  which  appear  to  be  the  only 
damped  media  for  which  energy  velocity  has  been  calculated  explicitly.®’ Lorentz 
media  are  defined^®  by  the  Maxwell  equations,  VxE  =  —dtB  and  VxH  =  dtD, 
and  by  the  constitutive  relations  B  =  fioH  and 

=  E{t)  +  f  d3g{s)E{t  -  3 

Jo 

g(t'j  =  sini/t. 

The  media  are  parameterized  by  a  plasma  frequency  o/p,  a  damping  constant  d, 
and  a  naturrd  frequency  t/.  In  the  frequency  domrun,  =  e{(jj)E{u)  and 

e(a;)  =  1  —  u>p/(a;^  —  w®  +  2iJa;)  (2) 

with  u;,  =  y/u^  +  iJ®.  Media  with  electric  dipoles  of  strength  that  oscillate  with 

damping  5  and  resonant  frequency  u/f  as  £  +  2Sx  +  lofx  =  constant  x  are 
Lorentz  media^^  in  the  sense  of  (1)  and  (2).  The  frequency  components  of  pulses 
in  Lorentz  media  travel  at  the  energy  velocity® 
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Ve  =  c  [(Re/:)  +  ^  k  =  Kek  +  ilmk  =  LJ^/e/c  (3) 

which  generalizes  group  velocity  in  the  sense  that^  a.s  5  0. 

Comparisons  of  group  and  energy  velocity  (Section  5)  may  be  the  parts  of 
this  paper  that  have  the  most  general  interest.  Figures  1  and  2  illustrate  one  such 
comparison  of  Ug  and  v^,  and  the  comparison  also  relies  on  Ref.  5.  That  paper 
computes  precursors  for  a  Lorentz  medium  (2)  whose  parameters  (wp,  Wr,  are 
on  the  order  of  10^®/sec.  For  one  computed  precursor,  the  instantaneous  frequency 
Winst  WcLS  measured  at  60  different  instants.  (Notice  that  measuTement  here  refers 
to  a  numerical  simulation,  not  a  laboratory  experiment.)  The  measurements  are 
graphed  in  Fig.  6  of  Ref.  5,  and  some  of  that  graph’s  features  are  sketched  here. 

All  60  mezisurements  lie  in  the  outer  dashed  regions  of  Figs.  1  and  2. 

Figure  1  shows  that  for  the  medium  (wp,  Wr,  J )  of  Ref.  5  and  Vg  for 
a  medium  (wp,  a;,,  0) — which  has  the  same  ujp  and  a;,,  but  has  5  =  0 — differ 
considerably  over  only  5-10%  of  the  frequency  range  for  which  Ref.  5  hais  measurements 
but  no  frequency  meaisurement  actually  lies  in  that  narrow  (5-10%)  range.  In 
the  complementary  region,  where  Vg  and  are  approximately  equal.  Fig.  6  of 
Ref.  5  has  six  unusual  measurements  that  seem,  at  first,  to  be  inconsistent  with 
both  Ve  and  v^]  but  those  six  frequency  measurements  are  unreliable’  for  reasons 
*  The  six  measurements  that  are  unreliable  for  our  purposes  are  useful  in  the 
context  of  Ref.  5.  That  paper  says  correctly  that  the  zigzags,  or  oscillations,  in 
the  graph  of  the  six  measurements  are  caused  by  interference  of  Sommerfeld  (high 
frequency)  and  Brillouin  (low  frequency)  parts  of  a  precursor.  The  paper  then 
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Fig.  1.  Precursor  measurements  in  Ref.  5  are  equally 
well  described  by  for  (u>p,  u>t>  S)  and  Ug  for  (wp,  Wr,  0). 
Group  and  energy  velocity  differ  significantly  near  the 
stop  band  (4xl0^®/sec,  6xl0^®/sec),  but  instantaneous 
frequency  was  not  measured  there. 
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Fig.  2.  Detail  from  Fig.  1.  Group  and  energy  velocity 
overlap  here,  to  within  the  width  of  a  pen  stroke.  In 
contrast,  many  measurements  in  Fig.  6  of  Ref.  5  deviate 
from  theoretical  curves  by  the  width  of  a  few  pen 
strokes,  owing  to  small  measurement  errors. 
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described  here  in  Section  5,  in  a  discussion  of  “zigzags”  and  multiple-valued  speeds. 
The  other  54  measurements  in  Fig.  6  of  Ref.  5  are  reliable  and  are  equally  well 
described  by  Vg  and  v^.  Figure  2  shows  the  high  degree  of  accuracy  with  which 
Ve  for  (ufp,  u>r,  S)  and  Vg  for  {a;p,  a;,,  0)  agree  over  the  (dashed)  range  of  reliable 
measurements.  This  paragraph  and  the  one  previous  explain  our  earlier  statement 

that  waves  computed  for  one  highly  absorptive  medium  are  described  with  equal 

« 

accuracy  by  group  and  energy  velocity.  Section  5  compares  Vg  and  further. 

The  oscillator  that  underlies  the  Lorentz  medium  studied  in  Ref.  5  has  a 
natural  frequency  u  that  differs  from  the  resonant  frequency  uij  by  0.25%.  That 
small  difference  accounts  for  the  similar  values  of  Vg  and  in  Fig.  2.  Increasing 
the  damping  S  would  broaden  the  frequency  range  over  which  for  (wp,  a;,,  E) 
and  Vg  for  (wp,  a;,,  0)  differ  significantly;  but,  in  a  competing  effect,  a  larger  damping 
would  also  severely  attenuate  a  broader  range  of  frequency  components.  With 
those  competing  effects  in  mind,  we  searched  for  Lorentz  media  in  which  Vg  and 
Ve  wotild  differ  significantly  for  some  frequency  components  that  could  be  observed 
and  mecisured.  Our  best  results — for  which  the  differences  between  Vg  and  v*  eire 
the  greatest  we  have  seen — are  in  Section  3.  The  section  describes  what  appears 
to  be  the  first  evidence  that  damping  affects  speed  noticeably,  though  the  effect 
is  slight.  Section  4  identifies  a  similar  small  effect  in  a  paper®  on  Debye  media. 
Section  5  reviews  comparisons  of  group  and  energy  velocity.  Section  6  introduces 
and  examines  an  elementary  method  for  measuring  permittivity.  The  elementary 
method  may  be  useful  when  damping  affects  wave  speed  slightly,  as  happens  in  all 
numerical  simulations  that  are  reported  or  cited  here, 
relates  the  onset  of  interference  to  a  prediction. 
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2.  INSTANTANEOUS  FREQUENCY 

Pleshko^^  was  the  first  to  measure  in  a  laboratory  the  instantaneous  frequency 
of  a  precursor,  which  is  a  frequency-modulated  wave.  Pleshko  assumed,  for  the 
salce  of  measurement,  that  consecutive  zero  crossings  of  a  time  trace  are  separated 
by  a  half-period,  as  is  true  for  a  sine  wave  of  fixed  frequency.  The  corresponding 
approximation  to  instantaneous  frequency 

f 

Winst  :=  7r/(r  -  0  (4) 

is  aissigned  to  the  instant  half-way  between  zero  crossings  T  >  t.  The  zero-crossing 
technique  is  not  always  useful.  For  example.  Figs.  3  and  4  have  no  zero  crossing 
after  2.3  x  10”^^ sec,  although  frequencies  there  obviously  decrease  at  those  late 
times. 

This  section  examines  measurement  techniques  in  which  extrema  and  inflection 
points  replace  zero  crossings  in  (4).  The  techniques  are  examined  carefully  because 
they  are  used  to  get  the  paper’s  main  result  (Section  3),  which  is  about  a  small 
effect.  Our  technique  for  measuring  instantaneous  frequency  consists  of  two  guidelines: 

(a)  Infer  frequencies  from  pairs  of  inflection  points 

or  pairs  of  extrema",  do  not  mix  types.  (5) 

(b)  Inflection  points  yield  better  accuracy. 

Use  of  technique  (5),  instead  of  the  zero-crossing  technique,  could  have  improved 
measurements  in  Refs.  5  and  7. 
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Fig.  3.  Time  trace  of  a  transmitted  field  produced  by 
a  sharp  incident  pulse.  The  medium  is  undamped. 
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This  section  tests  candidate  techniques  on  a  simulated  precursor.  The 
simulation  uses  computer  programs  from  Krueger  and  Winther  (personal  communication) 
that  are  described  in  Appendix  A.  The  simulation  is  for  a  3-mm-thick  homogeneous 
slab  of  Lorentz  medium  (1)  that  has  plasma  frequency  =  10^^/sec,  resonant 
frequency  a;,  =  5  x  10" /sec,  and  damping  J  =  0.  Outside  the  slab  is  free  space 

with  unit  permittivity.  The  incident  field  is  specified  just  inside  the  slab’s  left 

« 

WciU,  and  the  response  is  computed  just  inside  the  slab’s  right  wall. 

An  incident  delta  pulse  produces  a  transmitted  field  that  is  a  delta  function 
plus  a  remainder,  which  is  called  a  precursor.  The  precursor  is  illustrated  in  Fig.  3 
and  is  defined  analytically  in  the  paragraph  after  (A6).  We  use  inflection  points 
and  relative  extrema,  but  not  the  absolute  extremum,  to  measure  instantaneous 
frequency  in  three  ways.  Figure  4  focuses  on  the  useful  part  of  Fig.  3;  its  circles 
mark  inflection  points  and  its  asterisks  mark  extrema.  Inflection  points  in  Fig.  4 
do  not  foUow  the  slow,  upward  rise  expected  of  a  low-frequency  component;  yet 
we  win  see  that  inflection  points  give  the  best  reconstructions. 

Figure  5  illustrates  the  slowness  (c/v)  mecisured  from  Fig.  4  using  three 
techniques.  Circles  in  Fig.  5  come  from  substituting  consecutive  inflection  points 
for  consecutive  zero  crossings  in  (4);  the  asterisks  come  from  consecutive  extrema; 
and  diamonds  come  from  pairs  consisting  of  an  inflection  point  and  an  adjacent 
extremum.  The  curve  is  the  group- velocity  prediction  (c/ug)  for  the  medium, 
which  is  undamped.  Figure  5  is  numerical  evidence  that  the  inflection-point  technique 
(circles)  is  the  most  accurate  of  the  three  considered.  We  will  now  explain  this 
numerical  evidence. 
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Fig.  5.  Slowness,  as  inferred  from  Fig.  4  by  three 
different  methods.  The  curve  is  c/v^.  The  medium 
is  undamped. 
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The  signzJ  in  Fig.  4  has  a  high-frequency  oscillation  and  a  low-frequency, 
upward  drift.  A  low-frequency  drift  whose  effect  is  constant  between  two  extrema 
should  cause  little  error  in  measuring  u;in3t  because  extrema  are  zeros  of  the  first 
derivative,  and  the  derivative  of  a  constant  is  zero.  The  error  should  be  larger  for 
a  drift  that  is  linear  but  non-constant  between  extrema,  though  such  a  linear  drift 
between  inflection  points  (zeros  of  the  second  derivative)  produces  little  error. 
Thus,  frequency  mezisurements  based  on  inflection  points  should  be  more  accurate 
than  measurements  bcised  on  extrema,  and  far  more  accurate  than  zero-crossing 
measurements.  Further,  systematic  errors  from  extrema  and  from  inflection  points 
shoiJd  differ.  The  source  of  those  differing  systematic  errors  may  be  seen  in  Fig.  4, 
in  which  each  inflection  point  is  closer  to  the  extremum  on  its  immediate  left  than 
to  its  other  neighbor.  However,  some  of  that  shift  is  e.xpected  because  wavelengths 
increase  with  time. 

Instantaneous  frequency  is  the  reciprocal  of  a  difference  (4).  Therefore,  the 
effect  of  differing  systematic  errors  is  reduced  considerably  by  pairing  inflection 
points  with  other  inflection  points,  and  extrema  with  extrenia.  These  leist  three 
paragraphs  are  summarized  in  (5). 


3.  WAVE  SPEED 

This  section  examines  three  precursors  in  damped  media.  The  precursors’  wave 
speeds  are  measured  in  accordance  with  the  guidelines  in  (5).  This  section’s  last 
three  paragraphs  describe  the  paper’s  main  results. 
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Figure  6  shows  a  precursor  computed,  by  the  method  in  Appendix  A,  for  a 
half-space  of  Lorentz  medium  that  has  plasma  frequency  Up  =  10^^/sec,  resonant 
frequency  =  5  X  lO^'^/sec,  and  damping  J  =  2  x  10^  The  natural  frequency 

u  =  differs  from  the  resonant  frequency  by  8%.  The  incident  field  is 

a  square-modulated  sine  wave  that,  just  inside  the  material,  has  amplitude  10, 
has  carrier  frequency  o^c  =  9  X  10^^ /sec,  consists  of  ten  periods,  and  is  directed 
into  the  medium.  The  field  is  computed  for  a  depth  3  mm  inside  the  medium. 

In  the  simulation,  each  3-mm  spatial  interval  is  represented  by  700  points  and 
each  6-mm/c  time  interval  is  also  represented  by  700  points,  denning  a  space-time 
grid.*  The  simulation  was  compared  with  one  for  which  each  3-mm-by-6-mm/ c 
patch  is  represented  by  a  450-by-450-point  grid.  The  simulations  agree  until  1.6  x 
sec  but  disagree  thereafter,  owing  to  accumulated  discretization  errors.  In 
particular,  the  extrema  and  inflection  points  marked  in  Fig.  6  occur  at  the  same 
times  in  the  two  simulations,  to  within  the  resolution  of  the  coarser  (450  by  450) 
grid.  All  other  simulations  in  this  section  passed  similar  convergence  tests. 

Fxtrema  in  Fig.  6  are  marked  by  asterisks  ajid  inflection  points  by  circles. 
Instantajieous  frequency  [(4)  and  (5)]  and  arrival  time  axe  inferred  from  pairs 
of  inflection  points  and  from  extrema.  Each  frequency  component  is  assigned  a 
speed  t;  ==  3  mm/Tk,  where  is  arrival  time  at  the  3-mm  depth.  The  slowness 
c/v  of  each  component  is  compared  in  Fig.  7  with  three  slownesses  c/ve,  which 
correspond  to  energy  velocity  (3)  for  media  that  have  zero  damping  (solid  curve), 
*Even  though  At  =  2Ax/c,  the  simulation’s  evenly-staggered  grid  gij  —  (iAx,  (i-f 
2j)Ax/c)  has  slope  c  between  adjacent  points,  as  can  be  shown  by  graphing  the 
points. 
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•  DUT 


Fig.  6.  Time  trace  of  a  square-modulated  sine  wave 
after  transmission  through  a  damped  medium.  Numerical 
convergence  is  evident  until  the  last  circle  only. 
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the  correct  damping  (J  =  2  x  10^^/seCj  dashed  curve),  and  critical  damping 
(J  =  =:  5  X  10^^/sec,  dash-dotted  curve).  Data  from  inflection  points  are 

plotted  as  circles  and  those  from  extrema  as  asterisks.  The  effect  of  discretization 
error  is  illustrated  in  Fig.  7  wherever  an  error  bar  exceeds  a  point  marker’s  width. 

The  smallest  un-illustrated  error  bar  is  half  the  width  of  a  marker.  Vertical  error 
bajs  are  minute.  Appendix  B  defines  the  error  bars  and  Appendix  C  has  further 
error  analysis. 

The  most  error-prone  points  in  Fig.  7  are  not  useful  for  distinguishing 
among  energy  velocities.  Figure  8  focuses  on  the  useful  data.  The  slowness  inferred 
from  inflection  points  (circles)  is  consistent  with  energy  velocity  for  damping 
2  X  10^^ /sec  and  is  less  consistent  with  zero  damping.  Slowness  inferred  from 
extrema  (asterisks)  is  less  conclusive;  however,  we  showed  in  Section  2  that  data 
from  extrema  are  less  accurate  than  those  from  inflection  points. 

We  computed  a  second  precursor  for  the  medium  that  underlies  Figs.  6— 

8.  This  time  the  incident  pulse  is  a  delta  function.  (The  paragraph  after  (A6) 
discusses  delta  functions.)  The  response  3  mm  inside  the  dispersive  half-space  is 
a  delta  function  plus  the  remainder  illustrated,  in  part,  in  Fig.  9.  The  positions  of 
extrema  and  inflection  points  change  negligibly  when  the  number  of  grid  points  is 
reduced  by  60%,  implying  numerical  convergence  throughout  Fig.  9.  The  simulation’s 
absolute  minimum  (—50  units,  which  is  attained  at  10”^^ sec)  is  not  graphed  because 
it  is  not  used  to  compute  wave  speed. 

Slowness  for  Fig.  9  is  computed  with  inflection  points  and  extrema.  The 
result  is  Fig.  10,  in  which  circles  represent  data  from  inflection  points,  asterisks 
represent  extrema,  and  curves  represent  energy  velocity  for  the  same  damping 


II  -  15 


S I ouness  (c/v  ) 


F r  equenc  y  ( r  ad  / s  ec  ) 


Fig.  7.  Slowness,  as  inferred  from  Fig.  6,  for  which  S  = 
2  X  10^^/sec.  Measurements  are  consistent  to  varying 
degrees  with  c/ve  curves  for  the  underlying  .medium 
(dashed),  and  for  undamped  (solid)  and  critically 
damped  (dash-dotted)  media. 
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Slowness  (c/v) 


Fig.  8.  Detail  from  Fig.  7. 
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Fig.  9.  Part  of  the  time  trace  of  a  sharp  pulse  after 
transmission  through  a  medium  with  J  =  2  x  10^^ /sec. 
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Fig.  10.  Slowness,  as  inferred  from  Fig.  9,  for  which 
J  =  2  X  10^^ /sec.  Mezisurements  are  consistent  to 
varying  degrees  with  c/ue  curves  for  three  different 
media. 
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constants  as  in  Figs.  7  and  8.  All  figures  omit  error  bars  that  are  narrower  than 
point  markers. 

In  Fig.  10,  inflection-point  data  (circles)  are  consistent  with  energy  velocity 
for  2  X  10^^ /sec  damping,  except  possibly  for  the  slowest  component,  which  comes 
from  oscillations  in  Fig.  9  that  have  the  smallest  2implitudes.  Data  from  extrema, 
which  are  less  accurate  than  data  from  inflection  points,  zigzag  across  Fig.  10. 

Zigzags  are  most  likely  caused  by  inaccuracy  of  the  extremum  technique  (5).  (The 
scattered  asterisks  would  otherwise  imply  nonphysiczdly  that  some  frequency 
components  arrive  twice.  For  example,  the  component  w  =  1.3  x  10^^/sec  would 
have  speeds  v  =  c/1.75  and  v  =  c/1.9.)  Conversely,  the  smoothness  of  inflection- 
point  data,  with  a  single  exception  owing  to  oscillations  of  small  amplitude,  is 
further  evidence  of  those  data’s  reliability. 

The  last  simulation  uses  the  same  plasma  frequency  Wp  and  resonzint  frequency 
a;,  cis  before.  The  damping  S  is  increased  to  3  x  10^^/sec,  so  the  natural  frequency 
1/  dlifers  from  the  resonant  frequency  a;,  by  20%.  The  incident  ptalse  is  a  delta 
function  and  the  response  3  mm  inside  the  dispersive  half-space  is  a  delta  function 
plus  the  remainder  in  Fig.  11.  There  is  numerical  convergence  throughout  the 
simxilation.  Slowness  is  illustrated  in  Fig.  12,  in  which  the  dotted  curve  represents 
energy  velocity  for  3  X  10“ /sec  damping  and  other  notation  is  as  in  previous 
figures. 

We  have  now  compared  wave  speeds  of  three  precursors  with  predictions 
from  energy  velocity  (dashed  and  dotted  curves)  and  group  velocity  (solid  curves). 

No  single  data  set  is  dramatically  more  consistent  with  (with  damping)  than 
with  Vg  (no  damping.)  However,  the  evidence — Figs.  8,  10,  and  12 — is  more  conclusive 
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Fig.  11.  Part  of  the  time  trace  of  a  sharp  pulse  after 
transmission  through  a  medium  with  J  =  3  x  10“/sec, 
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Fig.  12.  Slowness,  cis  inferred  from  Fig.  11,  for  which 
5  =  Z  X  10^^/sec.  Measurements  are  consistent  to 
varying  degrees  with  c/u*  curves  for  three  different 
media. 


if  viewed  all  together,  and  is  strengthened  by  properly  weighting  the  data  from 
inflection  points  (circles)  and  extrema  (asterisks).  In  particular,  we  showed  In 
Section  2  that  asterisk  data  should  have  greater  error  than  circle  data:  a  conclusion 
supported  further  by  the  zigzagging  of  asterisks  across  Figs.  10  and  12.  Therefore, 
asterisk  data  should  have  less  weight  than  circle  data.  With  that  understanding, 
the  data  are  conclusive. 

There  is  convincing  evidence,  in  the  case  Up  =  10^^/sec  and  tJp  =  5  x 
10^^ /sec,  that  damping  speeds-up  frequency  components  noticeably,  to  a  degree 
that  is  consistent  with  energy  velocity.  The  effect  is  too  subtle,  or  our  measurements 
too  crude,  to  resolve  2  x  lO^Vsec  damping  from  3  x  10“/sec  damping,  though  the 
data  are  clearly  inconsistent  with  critical  damping  5  x  10^^ /sec.  That  is,  the  data 
are  in  rough  agreement  with  energy  velocity,  but  not  in  fine  agreement. 

The  numerical  data  have  three  sources  of  error:  discretization  error  in 
computing  precursors,  discretization  error  in  measuring  instantaneous  frequency 
(4),  and  crudeness  of  the  inflection-point  technique  (5).  The  first  two  sources  of 
error  are  analyzed  in  Appendix  B,  and  are  represented  quantitatively  by  error 
bars  in  several  figures.  Appendix  C  estimates  the  error  caused  by  the  inflection- 
point  technique  (5),  with  the  result  that  this  section’s  conclusions  are  unaltered. 

4..  DEBYE  MEDIA 

Albanese,  Penn,  and  Medina^  computed  precursors  for  water.  Their  model  is  a 
conducting  Debye  medium 
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cr  =  constant  e(w)  =  Ci  + 


a 


1  +  iuT 


(6) 


for  frequencies  from  zero  to  100  GHz.  (Reference  9  uses  ei  =  5.5£oj  o-  =  72.7, 

T  =  8.1  X  10~^^sec,  and  <r  =  10“®.)  Electrical  properties  at  frequencies  above  100 
GHz  are  modeled  differently.  However,  the  incident  field  is  a  square-modulated  1- 
GHz  sine  wave  of  ten# periods’  duration;  thus,  the  precursor  is  determined  almost 
entirely  by  the  Debye  model  (6).  One  of  the  precursor’s  most  striking  features 
is  that  it  hats  few  oscillations,  which  imply  that  wave  speeds  can  be  measured 
for  only  a  few  frequency  components,  and  those  speeds  are  approximately  equal. 

That  is,  the  damping  constaint  T  only  slightly  affects  the  speed  of  frequency  components 
that  are  observable  and  measurable  in  that  simulation.  Similar  effects — which  axe 
also  slight — are  seen  in  Section  3  and  Ref.  5. 

The  earlier  characterization  of  T  as  a  damping  constant  should  be  expladned. 

The  dzimping  constant  5  of  the  Lorentz  model  (2)  is  aptly  named  because  of  its 
role  in  the  model  equation  x  -t-  2Sx  -t-  =  constant  x  It  is  sensible  to 

call  T  a  damping  constant  for  the  Debye  model  because  the  permittivity  (6)  is 
a  small-o;  approximation^®  to  the  Lorentz  model  (2),  with  damping  T  =  25. 

Although  Lorentz  and  Debye  models  have  easily-identified  damping  constants, 
it  may  not  be  possible  to  extend  the  idea  of  a  damping  constant  to  all  dispersive 
media.  However,  there  is  evidence  from  three  independent  sources — Section  3  and 
Refs.  5  and  9— that  damping  affects  speed  slightly  or  negligibly  for  oscillator-type 
models. 
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5.  COMPARISONS  OF  GROUP  AND  ENERGY  VELOCITY 

A.  Open  Questions 

Energy  velocity  is  a  generalization  of  group  velocity,^  and  in  that  sense  Vc  is  better 
than  Ug.  That  is,  is  better  in  theory.  But  how  much  better  is  empirically? 

In  particular,  how  much  better  does  Ue  predict  the  speeds  of  frequency  components 
that  can  be  observed  and  measured? 

B.  Empirical  Similarities 

We  explained  after  (3)  how  and  Vc  predict  equally  well  the  speeds  of  frequency 
components  observed  and  measured  in  Ref.  5.  That  may  seem  puzzling  at  first, 
because  is  complex- valued  in  the  stop  bands^*^^  of  undamped  Lorentz  media, 
and  it  is  also  complex-valued  for  damped  media.  This  subsection  discusses  stop 
bands  at  length,  and  then  explains  how  and  can  sometimes  be  equally  useful 
when  there  is  damping. 

Energy  velocity  (3)  predicts  the  speed  of  all  frequency  components,  including 

those  in  stop  bcinds,  but  cannot  predict  speeds  in  the  stop  bands^‘^^  (a;^,  +  ^p) 

of  undamped  Lorentz  media  (a;p,  a;r,  0)  because  is  purely  imaginary  there, 

We  will  show,  however,  that  there  has  been  no  reliable  frequency  measurement 
in  a  stop  band;  consequently,  the  uselessness  of  Vg  in  stop  bands  has  no  proven 
empirical  consequence. 

It  appears  that  the  only  papers  that  measure  instantaneous  frequency  for 
Lorentz  media  are  the  present  one,  and  Refs.  5  and  7.  Figure  6  of  Ref.  5  shows 
that  the  precursor  there  has  no  observable  and  measurable  component  in  the 
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stop  band  (4  x  10^®/sec,  6  X  10^®/sec).  Components  in  the  stop  band  also  aren’t 
observed  here.  Thus,  only  one  set  of  stop-band  measurements  has  ever  been  reported. 
Those  measurements  are  in  Figs.  8,  10,  12,  and  14  of  Ref.  7  (called  [7,  Figs.]),  for 
which  the  stop  band  is  also  (4  x  10^®/sec,  6  X  10^®/sec);  The  subject  of  that 
paper  (signal  velocity)  and  this  one  (ug  and  v^)  differ.  However,  the  parts  of  [7, 

Figs.]  to  the  left  of  the  intervals  identified  there  with  the  mark  “5c >”  may  be 
relevant  here.  That  is  because  the  left-of-5c  parts  are  said,  in  text  that  follows 
equation  (2.12)  of  Ref.  7,  to  be  dominated  by  precursors.  Yet  we  will  see  that  the 
left-of-5c  stop-band  measurements  of  [7,  Figs.]  are  not  accurate  enough  for  our 
purpose  because  the  measured  speeds  are  multiple- valued  functions  of  frequency, 
and  because  of  other  matters  relating  to  energy  velocity.  This  is  explained  in  the 
next  two  paragraphs. 

Energy  velocity  (3)  is  a  single-valued  function  of  frequency.  That  is  why 
Section  3  says  that  the  zigzagging  asterisks  in  our  Fig.  10  represent  measurements 
that  are  too  inaccurate  to  be  useful  here.  For  example,  if  the  asterisk  measurements 
in  Fig.  10  were  reliable  (they  aren’t,  which  is  why  we  developed '(5)),  then  the 
component  a>  =  1.3  X  10^^/sec  would  have  speeds  v  =  c/1.75  and  v  =  c/1.9. 

In  another  example,  the  different  degrees  of  zigzagging  in  three  data  sets  in  the 
present  paper’s  Fig.  5  make  the  relation  of  multiple-valued  speeds  and  inaccuracy 
even  more  clear,  especiadly  in  light  of  the  text  in  Section  2.  Figure  5  suggests 
further  that  the  bandwidth  of  zigzags  is  comparable  to  the  error  in  measuring 
instzintaneous  frequency.  Thus,  the  degree  of  zigzagging  in  (7,  Figs.]  suggests 
that  the  figures’  left-of-5c  parts  have  data  in  the  stop  band  merely  because  of 
inaccurate  frequency  measurements.  Such  a  result  is  unsurprising  because  simple 
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differentiation  (see  the  sixth  paragraph  of  Section  2)  shows  that  the  zero-crossing 
technique  of  Ref.  7  is  especially  apt  to  be  inaccurate  when  frequency  components 
interfere,  as  they  do  in  the  simulations  underlying  [7,  Figs.]  That  interference, 
and  other  sources  of  measurement  error,  are  discussed  in  Ref.  7  on  pages  1438 
and  1439. 

We  have  just  seen  that  the  left-of-^c  stop-band  measurements  In  [7,  Figs.] 
axe  inconsistent  with  the  energy-velocity  prediction  that  speed  is  single-valued, 
but  those  measurements  are  inconsistent  with  in  other  ways.  For  example, 
energy  velocity  (3)  predicts  for  the  medium  of  Ref.  7  that  frequency  components 
in  its  stop  band  (4  X  10^^ /sec,  6  x  10^® /sec)  travel  no  faster  than  one-eighth  the 
speed  of  light  in  vacuum,  but  [7,  Figs.]  are  all  inconsistent  with  that  prediction. 

To  make  a  more  thorough  comparison  with  Vo  can  plot  (3)  directly  on  [7, 
Figs.]  and  then  observe  that  every  left-of-0c  stop-band  component  travels  faster 
than  Uc  allows. 

We  conclude  that  the  relevant  stop-band  measurements  in  [7,  Figs.]  are 
too  inaccurate  for  our  purpose*  because  the  corresponding  speeds  are  multiple¬ 
valued,  and  all  of  those  values  are  too  high.  The  only  other  papers  that  measure 
instantaneous  frequency  appear  to  be  the  present  one  and  Ref.  5,  but  neither 
has  stop-band  measurements.  Therefore,  the  inability  of  to  make  stop-band 
predictions  has  no  proven  empirical  consequence. 

*  More  precisely,  no  left-of-5c  measurement  In  Ref.  7  is  accurate  enough  that  it 
can  be  relied  on  to  be  in  a  stop  band.  The  measurements  may  still  be  useful  in  the 
context  of  Ref.  7,  which  emphasizes  signal  velocity. 
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In  a  fined  matter  concerning  Ref.  7,  we  want  to  make  it  quite  clear  that  the 
paper’s  time  traces  are  not  in  question.  The  time  traces  are  unquestioned  because 
the  paper’s  departures  from  the  predictions  of  energy  velocity  are  most  likely 
caused  by  inaccurate  frequency  measurements.  Some  inaccuracies  are  discussed 
in  Ref.  7  on  pages  1438  and  1439,  and  others  are  caused  by  the  zero-crossing 

technique’s  special  sensitivity  to  interference.  Although  better  frequency-measurement 

« 

techniques  are  now  available  in  (5),  Appendix  B  shows  that  the  error  in  measuring 

frequency  still  predominates:  Frequency  measurements  are  delicate. 

Let  us  now  consider  damped  Lorentz  media,  for  which  group  velocity  is 

« 

complex- valued.  A  complex- valued  Vg  seems  useless^  but  if  the  damping  S  affects 
the  speed  of  observable  frequency  components  negligibly,  then  v^.  for  (cjp,  0;^^  S) 
and  Vg  for  the  related  undamped  medium  (a;p,  0;^,  0)  differ  negligibly  in  their 
predictions.  That  is  exactly  the  case  for  the  medium  studied  in  Ref.  5;  thus,  as 
is  explained  sifter  (3),  Vg  and  are  equally  useful  there.  In  a  related  matter, 

Section  4  finds  further  evidence  that  damping  attenuates  the  components  it  speeds 
up,  maJcing  them  difficult  to  observe.  Those  competing  effects  of  speed-up  and 
attenuation  explain  how  Vg  for  a  damped  medium  and  Vg  for  a  related  undamped 
medium  can  predict  equally  well  the  speeds  of  observable  frequency  components. 

C.  Empirical  Differences 

The  present  paper  appears  to  have  the  only  empiriccil  evidence  that  is  better 
than  Vg.  In  particular,  Section  3  identifies  media  for  which  is  noticeably  better 
than  Vg,  though  the  difference  is  slight.  Efforts  to  find  Lorentz  media  in  which  v^ 
is  far  superior  to  Vg  were  unsuccessful,  but  not  all  Lorentz  media  were  examined. 
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In  spite  of  numerical  evidence  from  several  sources,  it  is  still  unknown 
whether  energy  velocity  is  any  more  than  slightly  better  than  group  velocity  from 
an  empirical  point  of  view.  That  uncertainty  is  striking,  considering  that  many 
papers  mention  energy  velocity. 


6.  ELEMENTARY  INVERSION 

f 

This  section  makes  an  observation,  apparently  unnoticed  since  precursors  were 
first 

predicted^^*"^®  in  1914,  that  wave  speed  determines  permittivity  elementarily. 
Limitations  of  this  elementary  inversion  are  discussed.  The  conclusion  is  in  the 
first  two  paragraphs  after  (8). 

It  is  easy  to  measure  group  velocity  for  undamped  media.  Substitution  of 
Vg  =  dxjjfdk  into  (7)  verifies  that  permittivity  can  then  be  computed  from  as 

e(a;)  =  u) 

where  n(u;o)  is  the  index  of  refraction  at  a  starting  frequency  u;o.  If  group  velocity 
is  measured  for  high  frequencies,  as  has  been  done  here  and  in  Ref.  5,  then  one 
could  choose  a;o  to  be  the  highest  such  frequency  and  assume  n(ti;o)  ^  1.  After 
all,^^  n(a;)  ->  1  as  w  oo.  Thus,  it  is  easy  to  measure  the  high-frequency  part  of 
permittivity  for  undamped  media.  The  same  measurement  appears  to  determine 
with  fair  accuracy  all  but  the  damping  constant  for  damped  Lorentz  media.  The 
next  two  paragraphs  examine  difficulties  faced  by  elementary  inversion  for  an 


r  du’ 

u,„n(u„)+cy 
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Isolated  interval  of  frequencies  that  are  too  low  for  the  approximation  n(a;o)  ~  1 
to  be  valid. 

Identifying  a  starting  point  (wq,  n(a;o))  in  an  isolated  band  of  low  frequencies 
is  problematic,  especially,  because  laboratory  apparatus  that  produce  precursors 
may  not  be  suited  for  measuring  an  index  of  refraction.  In  some  applications, 
however,  one  may  have  prior  knowledge  of  starting-point  data  to  within  reasonable 
error.  For  high-freqliency  starting  points,  the  error  is  small  and  the  reconstruction 
good.  For  low-frequency  starting  points,  starting-point  error  is  likely  the  predominant 
error.  This  predominant  error  causes  the  mathematical  difficulty  that  is  considered 
next. 

Assume  in  this  paragraph  that  e  is  causal,  that  group  velocity  Vg  is  known 
accurately  for  every  frequency  in  an  interval  [a,  b] ,  and  that  n  is  known  for  no 
point  in  [a,  6],  Set  £a  to  be  the  right  side  of  (7)  after  the  undetermined  product 
uon{utQ)  is  replaced  with  an  arbitrary  parameter  a:  ,• 


ea(ti')  := 


-2 


r  dcj'  V 

cc  +  c  — r-p- 

Ju.0 


(8) 


A  causal  permittivity,  such  as  e,  can’t  have  a  second-order  pole  at^^  a;  =  0.  This 
section’s  last  paragraph  proves  that  Ca  has  a  second-order  pole  at  a;  =  0  unless 
a  =  a;on.(a;o).  Thus,  there  is  a  unique  causal  which  is  e.  In  inverse  scattering, 
uniqueness  is  often  helpful,  but  not  so  here:  Anything  but  the  precisely  correct — 
hence  unlikely — choice  of  a  yields  noncausal,  nonphysical  permittivity.  Thus, 
elementary  inversion  for  an  isolated  band  of  low  frequencies  has  a  predominant 
error  that  ruins  causality. 
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We  conclude  that  elementary  inversion  (7)  can  saiely  be  used  for  data 
covering  an  interval  that  includes  high  frequencies.  Inversion  then  is  valid  throughout 
the  interval,  even  if  the  interval  also  includes  low  frequencies.  Elementary  inversion 
may  not  be  useful  for  isolated  bands  of  low  frequencies. 

Time-domain  experiments^^^^^  on  dispersive  materials  were  performed  in 
the  1960s  and  1970s  and  the  data  analyzed  in  the  frequency  domain  with  methods 
more  complicated  tlian  (7).  Time-domain  experiments  continue, partly  in  response 
to  new  methods  of  data  analysis^^"^^  in  the  time  domain.  Elementary  inversion 
(7)  appears  to  be  the  simplest  frequency-domain  inversion  available  for  time- 
domain  data.  Elementary  inversion  is  useful  to  varying  degrees  in  different  media. 

It  is  useless  for  the  Debye  medium  simulated  in  Ref.  9,  useful  for  Lorentz  media 
simulated  in  this  paper,  and  more  useful  for  the  Lorentz  medium  simulated  in 
Ref.  5.  For  example,  permittivity  can  be  computed  by  simply  integrating,  as  in 
(7),  the  Vg  measurements  in  Fig.  5  or  the  measurements  in  Fig.  6  of  Ref.  5. 

This  paragraph  proves  a  result  used  earlier  in  this  section.  We  assume,  for 
the  reason  given  after  (8),  that  neither  e  nor  Cq  has  a  pole  of  order  two  at  a;  =  0; 
then  we  show  a  =  cx;on(a;o).  Let  A  =  a  —  a;on(u;o).  Then  (8),  Vg  =  (kjfdk,  k  = 
(jJTi/c^  and  e  =  imply  ^  "h  where  the  function  A  =  A(A  +  2nu;)ci;  ^  also 

does  not  have  a  pole  of  order  two  at  u  =  0.  The  equation  A  =  A(A  -f  2ruj)u 
implies 


2  A^  u^A'^ 


(9) 
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Neither  A  nor  c  has  a  pole  of  order  two  at  a;  =  0.  Thus,  (9)  implies  A  —  0, 
finishing  the  proof. 


7.  CONCLUSION 

Section  3  shows  that  damping  can  affect  the  speed  of  frequency  components  noticeably, 
though  the  effect  is  slight.  The  effect’s  small  size  is  surprising  because  it  implies 
that  energy  velocity  is  only  slightly  more  accurate  here  than  is  group  velocity. 

Section  5  notes  that  there  are  now  four  independent  sources  of  evidence  that 
damping’s  influence  on  wave  speed  is  slight  or  negligible  for  some  Lorentz  and 
Debye  media.  One  interpretation  is  that  the  effect  is  small  because  of  two  competing 
mechanisms.  Specifically,  damping  attenuates  severely  the  frequency  components 
whose  speeds  are  most  changed. 

One  practical  consequence  of  this  work  is  that  the  elementary  inversion 
technique  in  Section  5  may  be  useful  for  some  damped  media,  as  well  as  for  undamped 
media.  The  technique  uses  wave  speed  measurements  that  include  high-frequency 
data.  Two  numerical  examples  are  mentioned  between  (8)  and  (9). 


APPENDIX  A:  DIRECT-SCATTERING  ALGORITHM 
This  paper’s  simulations  all  use  computer  programs  from  Krueger  and  Winther 
(personal  communication).  The  algorithm  that  underlies  the  programs  has  not 
yet  been  published,  although  it  was  rederived  in  detail  by  one  of  us  and  is  further 
validated  in  Ref.  24.  The  aJgorithm  uses  a  wave-splitting  technique  to  calculate 
time-domain  Green’s  functions.  (Wave  splitting  is  described  simply  in  Ref.  25.) 


Similar  time-domain  Green’s  functions,  which  we  CcJl  Krueger-Ochs  Green’s  functions, 
have  been  used  for  isotropic^®  and  anisotropic^^  nondispersive  slabs  that  are  continuously 
stratified,  and  for  cylinders  whose  dispersion  is  of  type^®  -^(^>0  =  + 

d3g{3)E{r,t  -  3)]  whe  re  r  is  radial  distance.  These  functions  have  also  been 
used  in  transport  theory.^® 

In  response  to  a  referee’s  suggestion,  we  document  the  unpublished  Krueger- 

« 

Winther  algorithm  by  describing  a  slight  generalization.  We  state  results  only. 

Details  axe  deferred  to  Ref.  24  and  a  possible  future  paper  by  Krueger  ajid  Winther. 

We  will  solve  Maxwell  equations  for  pulses  normally  incident  on  a  finitely 
thick,  flat  slab  with  properties”  '^xE  =  —fiodtH,  Vxif  =  J  -f  dtD,  J{z,t)  = 
cr{z)E{Zyt)f  and 


D{z,t)  =  £0 


E{z,t)+  j  dsg[s)E{z,t  —  s) 
Jo 


(>11) 


where  /xq  and  to  are  constants  and  g{t  <  0)  =  E[z  >  0,t  <  0)  =  0.  The  slab’s 
boundaries  are  z  =  0  and  z  —  L.  Let  c  =  l/^iioeo  and  define  dimensionless 
variables 


x  =  zlL  T  =  ctlL  .  (yl2) 

in  which  the  slab’s  boundaries  are  x  =  0  and  x  =  l.-The  wave  equation  for  the 
electric  field  is 

\dl  -  dz  +  B{x)dr]u{z,T) -r  dl  f  ds'y{s)u{x,T  -  s)  =  0  (.43) 

■  ”  '  Jo 
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where  B{x)  =  —p,oLca-[z{x)]  =  —fioLca{Lx),  7(r)  =  —Lc  =  —Lc~^g{Lc~^T), 

and  tt(i,  r)  =  E[z{x),  f(r)]. 

Define  convenient  basis  functions 


u 


dsd~u(^x^  s) 


(A4) 


The  functions  u~  are  bcisis  functions  in  the  sense  that  u  =  -r  u~ .  Outside  the 
slab,  where  5  =  7  =  0,  we  have^®  =  u'^(z  —  ct)  and  u~  =  ^“(2  ct);  thus, 
the  time-domain  boundary  condition  of  incidence  can  be  expressed  conveniently 
in  terms  of  and  u~.  The  wave  equations  for  are 


2(9x  -b  5r)u'*'  =  -2(5i  -  5r)Ti  =  B(l)  [u"*" (l ,  t)  +  U  (x.t)] 

+  dr  ^57(5)  [u'^(x,T  -  s)  -bU-(2,T  -  s)]  , 

Jo 


(A5) 


which  are  coupled  integrodirferential  equations.  The  last  paragraph  defines  Krueger- 
Ochs  Green’s  functions  G-(x,t)  and  a  function  a(i)  such  that 


.u'*'(x,r)  =  a(x)/(r  —  i) -b  f  dsf[3)G'^{x,T  —  3), 

Jo 


u  (x,r)  = 


[  dsf[s)G  (x,T  —  s) 

Jo 


iA6) 


solve  (A5)  subject  to  the  boundary  condicions  u  ■  {0,t)  =  /(r)  and  u-(x  >  0,r'  < 
0)  =  u“(x  >  l,r)  =  0,  for  which  the  incident  field  is  /.  Thus,  computing  G~ 
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solves  the  scattering  problem.  The  logical  flow  of  the  solution  is  (C?“,  /)  -4  — > 

u  =  u'^  +‘n~  E{zyt)  =  u{zfL,  ct/L). 

Equation  (A6)  and  the  boundary  condition  ^“(x  >  1,t)  =  0  show  that 
is  the  transmitted  field.  For  that  reason,  G‘^(l,r)  is  called  the  transmitted 
response  to  an  incident  delta  function  /(s)  =  ^(•s),  after  an  attenuated  delta 

function  a(l)J(r  —  1)  is  subtracted  from  the  transmitted  field  in  (A6).  .A.lso, 

« 

the  term  ds/(3)G'^(l,T  —  s)  is  called  the  precursor,  for  any  incident  field 
/  for  which  the  previous  integral  exists.  The  response  and  the  precursor  are 
computed  without  having  to  represent  delta  functions  numerically.  Indeed,  delta 
functions  are  used  here  only  in  physical  interpretations,  not  in  derivations^®  or 
numerics. 

Let 


a(x)  =  e.xp 


(A7) 


Krueger-Ochs  Green’s  functions  G—  satisfy  the  coupled  integrodifferential  equations 


2(5,  +  aOG+(x,T)  =  -2{d,-dr)G-{x,r) 

=  B{x)  [G'^(x,t)  -I-  G“(x,r)]  a{x)dr-t{x,r  -  x) 

+  [  ds  [9,7(x,  s)]  [G~(x,t  -  s)  +  G  (x.r-s)] 
Jo 

subject  to  boundary  conditions 


(.43) 
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G-(s,0+)  =  0 

G-{x,r<x)  =  0 

G'^{x,x'^)  =  a{x)  [  dxB^{x 

Jo 

^  G”(x,i+)  =  a(x)B(x)/4 

where  x"*"  stands  for  the  limit  r  |  x.  The  functions  G~  are  continuous  except 
as  described  in  (A9).  The  method  of  Ref.  30  can  be  used  to  prove  that  there  is  a 
unique  pair  of  solutions  G~  to  (AS)— (A9).  The  integrodifferential  equations  are 
solved  numerically  with  the  method  of  characteristics.  That  method  is  applied  to 
similar  equations  in  Refs.  26  —  29. 

APPENDIX  B:  ERROR  BARS 

This  appendix  is  written  mainly  in  terms  of  extrema.  The  results  also  apply  to 
inflection  points,  because  inflection  points  are  extrema  of  a  first  derivative. 

The  error  bcirs  in  Figs.-  7,  8,  10,  and  12  represent  the  roots-mean-square 
of  random  errors  in  Winst-  There  are  two  sources  of  random  error.  They  are  the 
inexact  computation  of  time  traces  and  an  inexact  method  of  identifying  extrema. 
The  next  two  paragraphs  estimate  the  roots-mean-square  of  those  two  types  of 
error.  The  root-mean-square  (rms)  of  the  superposition  of  the  two  errors  is  then 
calculated.  That  calculation  determines  the  rms  error  of  each  e.xtremum’s  time 
coordinate.  Consecutive  extrema  T  and  t  define  an  instantaneous  frequency/ 
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:=  >r/(r  -  1) 


(•81) 


and  the  errors  in  T  and  t  thereby  affect  a/jnst.  The  rms  of  the  resultant  error  in 
Winst  is  calculated  in  (B2),  and  is  used  to  define  error  bars.  In  this  paper,  each 
pair  of  error  bars  represents  a  computed  frequency  ttiinst,  plus-or-minus  the  rms 
error  in  Winst*  Notice  that  the  bars  represent  random  errors,  not  systematic  errors. 
Systematic  errors  are  discussed  in  Section  2,  where  it  is  shown  that  the  inflection- 
point  technique  has  smaller  systematic  errors  than  the  extremum  technique.  Appendix  C 
shows  that  systematic  errors  insignificantly  affect  the  agreement  pf  inflection- 
point  data  and  energy  velocity. 

Let  us  now  estimate  the  random  error  caused  by  inexact  computation  of 
time  traces.  The  estimate  is  based  on  convergence  tests  described  in  the  second 
paragraph  of  Section  3.  Recall  that  time  traces  are  computed  on  two  grids  whose 
resolutions  differ  significantly,  and  the  traces’  extrema  are  found  to  coincide  within 
the  resolution  of  the  coarser  grid.  That  convergence  test  shows  that  the  coarser 
computation  attains  what  is  commonly  called  convergence;  thus,  one  expects 
the  finer  computation  also  to  converge.  Consequently,  one  expects  the  extrema 
of  the  finer  computation  and  the  extrema  of  a  hypothetical,  exactly-computed 
time  trace  to  coincide  within  the  grid  spacing  A  of  the  finer  grid.  We  therefore 
assume  that  errors  in  the  time  coordinates  of  the  finer  computation’s  extrema  are 
distributed  uniformly  between  —  A  and  A.  The  corresponding  distribution  Pi{€) 

of  errors  ;  is  1/(2A)  for  |£|  <  A,  and  0  otherwise;  and  the  rms  error  dePj 
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is  A/ y/S.  That  is,  the  inexact  computation  of  time  traces  causes  errors  in  the 
time  coordinates  of  extrema,  and  the  rms  of  those  errors  is  approximately  A/-v/3. 

Let  us  now  estimate  the  rms  of  the  error  in  identifying  extrema.  This  new 
error  will  be  defined  by  describing  the  way  extrema  are  identified  in  practice:  A 
time  trace  is  computed  on  a  grid  of  time-resolution  A,  is  subjected  to  convergence 

test,  and  is  found  to  pass  the  test.  The  trace  represents  an  electric  field 

« 

at  times  zA,  for  various  consecutive  integers  z*.  Each  relative  extremum  z’o  A  is 
found  to  be  sharp  in  the  sense  that,  first,  |EA(ioA)|  >  max{|£'^([z'o-l]A)|,  |EA([io+ 
1]A)|},  and,  second,  the  numerically-computed  time  trace  has  only  one  relative 
extremum  at  each  broad  peak  in  Fig.  6,  9,  or  11.  The  sharp  extremum  is  then 
said  to  occur  at  time  z'o  A.  Figure  13  shows  that  error  can  result  even  if  the  resolution- 
A  computation  is  exact.  For  that  reason,  we  assume  that  the  error  illustrated  in 
Fig.  13  is  uncorrelated  statistically  with  the  error  that  is  estimated  in  the  previous 
paragraph.  Because  of  the  lack  of  correlation,  it  is  not  a  restriction  to  assume, 
merely  for  this  paragraph’s  error  estimate  and  for  calculations  leading  to  (B2), 
that  the  resolution-A  computation  is  exact:  that  is,  to  assume  that  each  computed 
value  JE^{zA)  equals  the  exact  electric  field  E{t)  at  t  =  zA.  Existence  of  dtE{t) — 
which  follows  from  the  constitutive  relation  (1)  and  the  Maxwell  equation  for 
dtD — then  Implies  that.E(t)  js  almost  symmetric  in  a  small  neighborhood  of  any 
relative  extremum  T,  as  illustrated  in  Fig.  13.  That  symmetry  helps  determine 
the  maximum  possible  difference  between  z'o  A  and  T,  given  that  z’o  A  is  a  sharp 
extremum.  The  maximum  difference  is  illustrated  in  Fig.  13,  in  which  I’oA  is  T  — 

A/2,  minus  an  infinitesimal,  and  Ex  (to  A)  is  infinitesimally  larger  than  Ex  ([z’o  + 

1]A).  In  fact,  z'oA  is  a  sharp  e.xtremum  if,  and  only  if,  it  is  in  the  interval  (T  - 
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Fig.  13.  Extrema, (*)  and  grid  points  (•)  are  unlikely 
to  coincide. 
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A/2,T  +  A/2).  The  corresponding  distribution  P2(^)  of  errors  €  is,  consequently, 

I'/A  for  le]  <  A/2,  and  0  otherwise;  and  the  rms  error  is 

A/v/l2. 

The  roots-mean-square  of  two  sources  of  random  error  have  now  been  estimated. 
One  rms  error  is  Afy/Z  and  the  other  is  A/\/T2.  The  second  estimate  comes 
from  a  model  in  which  the  two  sources  of  error  are  assumed,  with  justification, 
to  be  uncorrelated.  Consequently,  the  rms  of  the  total  random  error  in  the  time 
coordinate  of  any  extremum  is  [(A/\/3)^  +  (A/a/T2)^)^^^  =  A-y/o/12  ^  2A/3. 

The  rms  error  that  was  just  calculated  affects  the  time  coordinates  T  and 
t  of  consecutive  extrema,  which  define  the  instantaneous  frequency  in  (Bl). 

A  relevant  equation  for  the  propagation  of  rms  errors  is  known  most  widely  in  the 
context  of  Gaussian-distributed  errors,^^  but  the  equation  also  is  valid  for  other 
distributions.^^  That^  propagation-of-error  equation  implies 

ii..ii=  (ii^Tii^)' +  (iie.ii^)'  '  m 

where  ||e^,||  is  the  rms  error  in  ix;in3t,  Ikrll  =  lkt!|  =  Ai/5/12  are  the  rms  errors  in 
T  and  t,  and  the  partial  derivatives  in  (B2)  are  calculated  by  differentiating  (Bl). 

In  this  paper,  error  bars  span  the  frequency  intervals  (u?;n3  t  ll^uilj  ^'inst  ~  ll^wlD* 

The  analysis  in  the  previous  paragraph  can  be  adapted,  with  changes,  to 
the  vertical  coordinates  cjv  =  c(T  -f  t)/(2' depth)  of  Figs.  7,  8,  10,  and  12.  The 
adapted  analysis  shows  that  vertical  error  bars  are  much  smaller  than  the  height 
of  the  figures’  point  markers.  Thus,  the  time  traces  in  Figs.  6,  9,  and  11  come 
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from  convergent  numerical  computations,  and  the  corresponding  errors  in  cfv  are 
minute,  yet  there  are  significant  errors  in  merisuring  frequency.  This  observation 
supports  the  statement,  made  in  the  next-to-last  paragraph  of  subsection  5B,  that 
frequency  measurements  are  delicate. 

APPENDIX  C:  P,ERTURBATION  OF  INFLECTION  POINTS 
Section  2  shows  that  the  inflection-point  technique  yields  our  most  reliable  data, 
but  those  data  still  have  errors.  Those  errors  are  estimated  here  in  two  ways.  The 
estimates  do  not  change  the  paper’s  conclusions. 

Inflection  points  follow  imperfectly  the  slow  upward  drift  in  Fig.  6.  There 
is  an  urge  to  displace  points,  especially  the  negative- valued  one  at  1.2  x  10“^^sec, 
to  slightly  later  times:  displaced  points  would  follow  the  drift.  To  estimate  the 
effect  of  such  displacement,  imagine  that  only  the  most  deviant  point — the  one  at 
1.2  X  10“^^sec — were  moved.  Then  the  circle  in  Fig.  8  at  2.2  x  10^^/sec  would 
draw  toward  a  slightly  lower  frequency  and  the  circle  at  1.6  x  10^^/sec  would 
move  to  a  higher  frequency,  improving  the  fit  to  2  x  10“ /sec  damping.  Thus, 
replacing  inflection  points  with  points  that  more  closely  follow  the  drift  in  Fig.  6 
would  reinforce  the  paper’s  conclusion  slightly. 

The  next-to-Iast  paragraph  of  Section  2  notes  a  second  suggestion  in  the 
data  that  inflection  points  be  moved  to  later  times.  In  the  context  of  that  paragraph, 
systematic  error  can  be  estimated  pessimistically  by  shifting  two  inflection  points 
in  Fig.  6  (the  ones  at  1.2  x  10““sec  and  1.4  x  10~“sec)  by  equal  amounts, 
to  about  the  midpoints  of  neighboring  extrema.  Such  a  shift  would  not  affect 
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instantaneous  frequency  (4)  but  it  would  increase  slowness,  which  is  1.3,  by  1%. 

A  1%  shift  of  the  corresponding  point  in  Fig.  8  would  not  bring  it  into  line  with 
zero  damping.  A  similar  shift  of  the  inflection  points  that  correspond  to  slowness 
1.5,  would  be  less  justifiable  because  a  neighboring  extremum  is  outside  the  range 
from  0  to  1.6  x  10“^^sec  for  which  convergence  in  Fig.  6  is  evident.  In  any  case, 
the  estimated  shiTc  for  slowness  1.5  would  not  bring  the  corresponding  point  in 
Fig.  8  into  line  with  zero  damping. 

We  have  just  shown  that  error  from  inflection  points  is  acceptably  small. 
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Inversion  Method  for  Depth-dependent  Dispersive  Media 


« 


1  Problem  Statement 

We  wish  to  use  reflection  data  to  determine  the  dispersive  properties  of  a  flat  slab  that  has 
boundaries  at  i  =  0  and  £.  The  slab’s  dispersive  properties  are  defined  by  the  Maxwell 
equations  V-D  =  V-5  =  Vxi?  +  d.fl  =  0,  B  =  fioH,  and  Vxif  =  dtD  +  J,  and  by  the 
constitutive  relations  [l]-[2] 

D(x,  t)  =  eoE{x,  t)+  f  djWb(s,  t  -  s) 

(1)  J. 

J(x,t)  =  ff(x)£(2,0+  /  d3Wi{x,t-3)Eix,s). 

That  two-sentence  problem  statement  assumes  that  the  permittivity  eo  and  permeability  po 
of  the  slab  are  equal  to  their  free-space  values,  as  is  true  for  many  nonmagnetic  materials. 
The  slab  (1)  is  stratified  in  one  spatial  dimension  x,  and  the  dispersive  properties  in  (1)  are 
functions  of  the  variable  t  also.  Such  functions  of  the  two  variables  x  and  t  will  be  inferred 
from  reflection  data  as  a  function  of  the  time  t  and  the  angle  of  incidence  3.  We  will  now 
use  a  gauge  for  dispersion  to  simplify  (1)  and  to  establish  a  fundamental  physical  limitation 
on  the  ability  of  purely  electromagnetic  experiments  to  determine  <r,  Wp,  and  Wj. 

A  dispersion  gauge  was  developed  recently  [3];  it  is  easy  to  understand.  A  gauge  is  a 
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free  parameter  that  does  not  affect  fields.  In  electrostatics,  for  instance,  the  notation 

(2)  ^(a)  =  <f>  +  a 

reminds  us  that  the  scalar  potential  is  arbitrary  to  within  an  additive  gauge  constant  a 
because  the  electric  field  E  =  -  V(^  +  a)  is  independent  of  a.  For  that  reason  the  potential 
is  said  to  be  gauge  dependent  and  E  is  said  to  be  gauge  independent.  It  is  physically 
impossible  to  measure  gauge-dependent  quantities:  That  is  why  two  potentials  and  i>(a) 
can’t  be  me^lsu^ed  separately,  although  the  difference  -  V'(a)  =  {<f>  +  a)  -  {i>  +  a)  =  <j>  —  r{i 
is  gauge-independent  and  physically  measurable.  Many  undergraduate  [4]  and  graduate 
[5]  physics  texts  use  the  vector  potential  as  a  second  example  of  the  physical  impossibility 
of  measuring  gauge-dependent  quantities.  This  physical  principle,  that  gauge-dependent 
quantities  can’t  be  measured  directly,  also  applies  to  the  dispersion  gauge. 

The  dispersion  gauge  is  any  function  f{x,t)  that  has  a  time  derivative  and  for  which 
/(i,  0+)  =  lini«;o  /(*,  t)  exists  and  is  finite.  The  gauge  affects  the  displacement  D,  the  current 
J,  and  the  dispersion-related  quantities  <t,  Wdi  and  Wj  by  way  of 

VxH  =  dtD[f)  + Jy) 

•0(/)(2iO=  £oE{x,t)+  I  d5WD,(/)(®,i- •»)JS;(2,5) 

Jo 

(3)  Ju){xA)  =  cry){x)E{x,t)  + 3)E{x,3) 

W"D.(/)(x,t)  =  Wb(x,i)-/(i,0 
W3,^f)ix,t)  =  W3(x,t)  +  dtf{x,t) 
er(y)(s)  =  ^(z)  + /(s,  0-1-). 

It  is  easy  to  verify  [3]  that  H  and  E  are  independent  of  /;  that  is,  H  and  E  are  independent 
of  gauge.  We  will  now  show  that  the  totaJ  current  J  and  the  dispersion  kernel  7 

(4)  J{x,t)zz  ^tD^J){x,t)-i■J^f){x,t) 

(5)  7(as,  t)  =  cr(y)(*)  -f  W’d,(/)(x,  t)+  f  daWj.(/)(x,  j) 

Jo 

are  also  independent  of  gauge.  The  total  current  J  is  independent  of  the  gauge  /  because 
the  gauge  independence  [3]  of  H  implies  that  both  sides  of  the  first  equation  in  (3)  are 
independent  of  /.  The  gauge  independence  of  7  is  established  by  using  (3)  to  expand  the 
terms  on  the  right-hand  side  of  (5).  We  have  just  established  a  fundamental  result: 

Theorem  1.1  Each  of  the  quantities  •D(/),  ^(/)»  WD,(y)>  und  is  gauge  depen¬ 

dent.  The  fields  E  and  H,  the  total  current  J,  and  the  dispersion  kernel  7  are  gauge  inde¬ 
pendent.  Consequently  7(1, 0-f)  and  dty{x,  i)  also  are  independent  of  gauge.  If  Wb  and  W3 
are  sufficiently  smooth,  theny{x,0+)  =  cr(x)  + Wd(x, 0+)  and  dty{x,t)  =  W3{x,t)  +  dtWx3{x,t). 

This  elementary  theorem  is  important  because,  as  the  second  paragraph  of  this  section 
reminds  us,  gauge-dependent  quantities  can’t  be  measured.  That  is,  in  the  framework  of 
the  constitutive  relations  (1),  it  is  impossible  to  measure  the  gauge-dependent  quantities 
cry),  WD,y),  and  T/j,(y)  separately;  for  that  reason,  the  following  sections  of  this  paper 
concentrate  on  measuring  the  gauge-independent  combination  dty  =  9tWD,(;)  + 

A  Cautionary  Note:  Theorem  1.1  (above)  has  the  following  cautionary  implication 
relating  to  energy  density.  Many  energy-density  computations  are  btised  on  Poynting’s 
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equation  E  •  (VxH)  =  E  •  J,  in  which  both  sides  are  independent  of  gauge.  However, 
some  matters  associated  with  the  related  energy-density  computations  are  “open  to  some 
criticism,”  says  [6]  in  a  discussion  of  nondispersive  media;  but  for  dispersive  media  there  is 
one  more  matter  of  concern:  Some  energy-density  computations  [6]  have  separate  physical 
interpretations  for  the  terms  E-dtD^j)  and  E  •  J[j)]  yet  theorem  1.1  (above)  implies  that 
those  terms  can  not  be  measured  separately.  The  term  E  •  J*,  however,  is  independent  of 
gauge,  and  can  therefore  be  expected  to  be  physically  measurable. 


2  Method  of  Solution 

Apply  the  gauge  /  =  Wb  to  (3)  to  obtain  the  dispersive-conductivity  representation 

J{x,t)  =  a{x)E{x,i)+  [  dsg{x,i-  3)E{x,3) 

(oj  Jo 

a(x)  =  7(2, 0+)  =  cr(x)  +  Wb(a!i  0+) 
g{x,i)  =  dt^{x,t)  =  Wj(x,t)  +  diWu{x,i) 

of  the  constitutive  relations  in  (1).  Theorem  1.1  says  that  5  and  g  are  gauge ‘independent. 
One  published  assertion  [7]  may  be  interpreted  as  suggesting  that  a  must  be  0;  but  the 
assertion,  as  printed  in  [7],  is  not  substantiated  in  the  text  of  [7].  In  this  presentation 
we  assume  merely  that  a{x)  is  a  known  function;  we  do  not  necessarily  assume  a  =  0. 
Our  inverse  problem,  therefore,  is  to  infer  g{x^t)  from  time-dependent  reflection  data  as  a 
function  of  the  angle  of  incidence  9. 

We  assume  that  5  €  Loo[0,  L]  and  g  G  ioo([0,  L]  x  [0, 00)).  Then  the  techniques  of  [8]  and 
[3]  imply  that  there  are  integrable  functions  R  and  T  such  that 


-Xf  ft^hjco 

{Tf){i)  =  «P  -  y  dx5(i)/(2co)  fit  -  L/co)  +  d:,T{i  -  L/co  -  s)f(s) 


are  the  time-dependent  reflected  field  and  transmitted  field  {Tf){t)  produced  by  a 

field  /  incident  on  the  spatially  heterogeneous  slab  (6)  of  length  L.  The  speed  of  light  in 
free  space  cq  =  Ijy/poeo  is  used  in  (7). 

The  Redheffer  operator  product  [9]  shows  that  the  reflection  operator  for  a  composite 
formed  of  two  slabs  is 


(8)  ■  7^c  =  7^l  +  rl(l-7^27^l)-^7^2Tl, 

where  each  numerically  subscripted  operator  has  the  form  (7).  A  little  thought  shows  that 
the  inverse  operator  in  (8)  exists  and  has  a  resolvent-kernel  representation  [lO]  that  is  related 
to  a  Volterra  integral  equation  of  the  second  kind,  whose  difference  kernel  is  *  ds‘R2{t  — 
5- j')Ri(a').  The  resolvent-kernel  representation  and  the  integral  representations  in  (7)  are 
used  to  expand  the  four-operator  composition  Ti(l  —  'R2R-\)~^R-2T\  in  (8)  in  terms  of  the 
integral  kernels  of  the  operators.  A  discretization  of  that  four-operator  composition  is  used 
to  solve  the  inverse  problem. 

A  long  derivation  shows  that  the  integral  kernel  Ri  of  the  operator  Tic  has  the  form 

(9)  Ri{t  -1-  ^)  =  A{{t)  +  +  <!)-}-  Ct{t)Rl{t  -l-  5  -  2xi/co)  -1-  0(J3), 
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in  which  /ic(t  +  ^)  is  the  reflection-kernel  data  at  a  time  t  +  5  that  is  slightly  later  than  t, 
Gi(t  +  J)  =  ff[aso.(t  +  ^)]  is  evaluated  at  the  first  spatial  discretization  point  xo  =  0,  Rl{t + 
J  -  2xi/co)  is  the  reflection  kernel  for  the  whole  slab  with  the  layer  (0,  ii)  removed,  and 
the  quantities  A{,  B(,  and  C®  depend  on  g  for  times  equal  to  or  earlier  than  t  and  on 
for  times  equal  to  or  earlier  than  t  -  2xi/co.  The  superscripts  in  (9)  are  reminders  that, 
although  scattering  is  generally  a  function  of  the  incident  angle  9,  the  material  properties 
g{x,t)  are  independent  of  9.  The  computation  of  A®,  Bf,  and  Cf  is  described  in  the  next 
section. 

The  next  paragraph  shows  how  (9)  is  used  for  inversion  when  the  slab  [0,  L)  is  discretized 
into  three  spatial  intervals.  The  three-interval  discretization  is  generalized  to  n  layers  later 
in  this  section. 

This  paragraph  describes  the  central  idea  in  the  new  inversion  method  used  here.  This 
central  idea  is  illustrated  by  discretizing  the  slab  [0,Zr)  into  three  spatial  layers  [0,  xi), 
[ii.zj),  and  [xj.I),  and  adding  a  free-space  layer  [i,oo).  Omitting  time  dependences  such 
“  Wi  (*  +  <^)i  (t -1-  <1  -  2i,7co),  and  so  forth,  we  iterate  (9)  as  follows: 


Ji|  =  A?  +  BfGi  +  C{(,A\  +  BiGi  +  ClR-i)  = 

(10)  it®  =  A®  +  B®Gi  +  G®A®  +  G®5®G2  +  G®C®(A®  +  B®G3  +  Cj®/!®)  =' 
iti  =  (A®  +  G®A®  +  G®G®A®)  +  B®Gi  +  GfSfGj  +  C®G®B®G3. 

The  recursive  derivation  (10)  follows  from  (9)  in  the  context  of  layer  stripping.  In  particular, 
(9)  relates  the  reflection  itf  from  the  four-layer  composite  (0,  co)  =  [0,  ii)  U  [xi,  xj)  U  [xj,  L)  U 
[2i,  oo)  to  the  dispersion  kernel  Gi  of  the  first  layer  [0,  xi)  and  to  the  reflection  it®  from 
the  three-layer,  stripped  composite  [xi.co)  =  [xi.xj)  U  (x2,i)  U  (ii,oo).  But  the  three-layer 
composite  can  be  similarly  related,  by  using  (9)  again,  to  the  reflection  from  the  two- 
layer,  stripped  composite  [12,00)  =  [xi.ii)  U  [i/,00).  Thus,  one' iteration  of  (9)  yields  the 
first  line  of  (10),  in  which  it®  is  the  reflection  kernel  for  the  two-layer,  stripped  composite 
[x2,L)  U  [ii,co).  The  recursive  derivation  (10)  terminates  with  the  reflection  kernel  itj  =  0 
because  the  reflection  it®  from  the  free-space  layer  [L,  00)  is  zero.  Because  of  the  recursive 
nature  of  the  derivation,  the  numerical  algorithm  that  is  used  to  compute  the  coefficients  A®, 
J5f,  and  Cf  for  the  first  composite  [0, 00)  can  also  be  used  to  compute  the  general  coefficients 
A®,  5®,  and  C®  for  any  of  the  stripped  composites  [x,_x,oo).  The  algorithm  for  computing 
A®,  B®,  and  Cf  is  described  in  the  next  section.  Please  note  that  every  5-dependent  term 
in  (10)  is  marked  with  a  5  so  as  to  emphasize  that  the  5-independent  unknowns  Gi,  G2, 
and  G3  can  be  determined  by  measuring  the  reflection  Ri  from  the  full  composite  [0, 00)  for 
three  angles  of  incidence  5,  and  then  inverting  the  three-by- three  matrix  equation  that  is 
implicit  in  the  last  line  of  (10). 

We  will  now  present  a  matrix  equation  for  an  n-layer,  (n+l)-point  spatial  discretization 
xo,  xi,  . . i„  of  the  slab  [0,  L).  The  n-layer  generalization  of  (10)  is 

(11)  B{Gx  +  B^.Gi  n  Ci  =  B®  -  A®  -  AJ  n  G®, 

>  =  2  Jk=l  i  =  2  k  =  l 

in  which  all  the  Gj  terms  are  on  the  left  side  of  the  equality.  If  is  measured  for  N  angles 
9  then  (11)  is  equivalent  to  the  matrix  equation 

(12)  M-[Gl,G2,...,Gnr  =  [t^1.V2,...,V;/P. 
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for  which  the  i,  j  element  of  the  matrix  M  is 


(13) 


= 


1  111=1 


ifj  =  1; 
otherwise 


and 

(14) 


Jk-l 


Vi 


—  Ri^c  ~  «^»,1  ^  ^  J][ 


k  =  2  1=1 


and  denotes  transposition.  For  each  i  =  the  quantities  Aij^  Bij,  Cij^  and 

Ri^c  stand  for  Cy,  and  evaluated  at  the  incident  angle  9  -  9i.  If  the  number 

of  measurement  angles  N  is  greater  than  the  number  of  layers  n,  then  the  system  (12)  is 
overdetermined.  If  iV'  =  n  then 

(15)  [Gu  . . . .  GnF  =  M- 1  .  [vi,  ,  VnP. 


3  Computational  Algorithm 

An  inverse  scattering  problem  for  heterogeneous  dispersive  media  is  defined  in  (6).  The 
last  steps  in  the  algorithmic  solution  are  (12)-~(15).  All  that  remains  is  to  show  how  the 
quantities  Aj,  and  Cj  in  (12)“*(15)  are  related  to  the  material  parameters  defined  in  (6). 
That  relation  is  presented  here  in  the  form  of  a  computational  algorithm.  The  algorithm 
comes  from  derivations  that  are  much  too  long  to  report  here. 

The  following  three-page  computational  algorithm  documents  the  extent  of  development 
of  a  solution  to  the  inverse  problem  identified  in  this  paper’s  first  paragraph.  The  three- 
page  algorithm  itself  contains  no  new  ideas;  the  new  ideas  in  this  paper  are  all  in  section  1 
and  in  the  paragraph  that  includes  equation  (10).  The  algorithm  is  worthwhile  because  it 
is  ready  to  be  programmed. 

First  define  c  =  [/iofo-(sm^  ^)/co]"^/^,  (  =  ac/i,  and  r  =  ci/i,  in  terms  of  the  total  length 
L  of  the  slab  and  the  angle  of  incidence  9,  The  slab  interval  is  [0,  L)  in  the  ordinary  spatial 
coordinate  x  and  [0, 1)  in  the  unitless  coordinate  In  this  section,  functions  of  unitless 
variables  such  as  5(^,  r)  are  understood  to  be  related  to  functions  of  x  and  t  according  to 
the  rule  5(^,r)  =  g{^L,TL/c).  We  will  discretize  the  unitless  variables  as  ^  =  i/i  and  r  =  2jh^ 
for  any  h  that  is  the  reciprocal  of  a  positive  integer. 

For  the  sake  of  computational  efficiency,  in  case  the  slab  has  some  thick  homogeneous 
layers,  we  define  indexes  0  ^  Lq  <  Li  <  L2  <  • '  -  <  Ln  —  I//1  for  the  boundaries  of  the  layers; 
in  particular,  for  each  i,  5(^,r)  and  5(0  are  constant  €  [hLi^i.hLi).  The  layer  indexes 
do  not  restrict  the  generality  of  this  work  because,  if  there  is  no  thick  homogeneous  layer, 
then  we  set  I*  =  i  for  i  =  0, 1,2,  ...,n.  The  notation  L{i)  =  U  is  used  later  in  subscripts 
such  as  gn  -  9 Hi)- 

The  algorithm  in  this  section  is  derived  using  layer  stripping,  as  is  symbolized  in  (10). 
If  we  strip  layers  1  through  (i  -  1)  from  the  laminated  slab,  then  the  remainder  1) 

of  the  laminate  is 


where  H  is  Heaviside’s  step  function.  Each  stripped  layer  [kLi^i.hLi)  has  the  electromag¬ 
netic  properties 

9\:hi,  r)  =  T)H{hLi  -  hLi.r  -  0 
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Let 

(18)  =  g{hLi.u  2jh),  5(0  =  cF(/ii._i) 

represent  the  slab’s  properties  at  the  left  face  of  the  i***  layer  [hLi-i,hLi),  and  let 

(19)  p(*)  =  — /ioA(Xj  —  X,-_i)£ca(‘)/2. 

Let  Gr^'’^(^,t)  be  Krueger-Ochs  Green  functions  [8],  such  that 


JO 


dsdiE^^\^,3) 


0 

(0-/ 


(20) 


4^'*’U.t)  =  exp  [-3?(OeL/(2c)] 


f{r-()+  r"d.Gij)+(^.r -.)/(,) 

VO 


/  d5G|(^"(^,r- j)/(3) 

Jo 


solve  the  scattering  problem  of  fields  /  incident  on  the  stripped  slab  (i  =  1, 2, ....  n). 

Similarly,  let  be  the  Krueger-Ochs  Green  functions  that  solve  the  scattering 

problem  for  the  isolated,  stripped  layer  (i  =  1,2,...,  n).  Such  Krueger-Ochs  Green 

functions  were  originally  developed  analytically  and  numerically  for  nondispersive  media 
[8],  then  generalized  and  numerically  validated  for  normal  incidence  on  dispersive  media 
[11],  and  then  generalized  to  oblique  incidence  [3]. 

Let  us  now  re-introduce  the  time  arguments  that  were  omitted  from  (10)-(15).  Equation 
(11)  then  becomes 

n  j-l 

B\g{hLc,T  +  <y)  +  X! r  +  S-  2hLj-i)  = 


(21) 

J=2 

n 

Jk  =  l 

i-1 

=  i2.^(r  +  5)>l?(r)-X>l,V- 

2hL,..x)]lCt 

k=l 

Define 

4‘>  =  At[2{j-U.,)h] 

SCO  =  Bf 

(22) 

C(‘)  =  Gf 

G^l'!t  =  GP^iih,2Jh), 

with  4’^  =  0  when  j  <  X<_i.  Replace  r  and  S  with  2Jh  and  2h  in  (21)  to  obtain 
(23)  BMg%  +  X  f[  G(')  =  Rt[2{J  +  l)h]  -  A^p  ~ X  11 

k  =  2  /=:!  ik=:2  /=! 

Equation  (23)  is  of  the  same  form  els  (11),  except  that  (23)  denotes  the  time  dependence 
more  explicitly.  Equations  (21)-(23)  indicate  correctly  that  and  Cl*)  are  time  indepen¬ 
dent;  in  fact,  a  long  derivation  shows  that 

s(‘)  =  (A/2)[2  +  Ap(‘)]/  [2  -  hpl')  - 


(24) 


C(*)  —  exp  [p(*)Aj/c]  I  1  +  [^L,ko*^R,0,0  ^L,A(i).0^L,A{i).0 
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where  A(;)  =  Aj  =  X,+i  —  Li. 

The  only  task  remaining  is  to  show  how  to  compute  That  task  is  the  longest  one 
in  this  section  and,  because  of  a  triple  integral  in  (27),  it  is  also  the  most  computationally 
intensive  task  here.  The  triple  integral  is  of  convolutional  form,  however,  so  the  fast-Fourier- 
transform  (fft)  technique  for  computing  convolutions  can  be  applied.  The  fft  technique 
would  speed  up  the  triple-integral  computation,  but,  if  one  of  the  integrands  has  a  Fourier 
transform  that  is  broad  banded,  then  the  fft  technique  may  be  significantly  less  accurate 
than  would  a  straightforward  trapezoid- rule  computation  of  the  triple  integral  in  (27). 

Define 


(25) 


KW{t)*=  fds  t  -  ,)  G<;)-  (0, 5) 

Jo 

p(^{t)  =  +  J  3)pi'){3) 

77(0(0  =  p(*)(t)  +  Gp(/iA.-,t) 


p(0 


exp[khpi')]gj 


(«•) 


j-i 


J  =  l 


The  quantity  p(0  is  the  resolvent  kernel  [lO]  for  a  Volterra  integral  equation  of  the  second 
kind,  whose  difference  kernel  is  3)  =  if  5).  A  time-discretized  version  of  can 

be  computed  by  inverting  a  triangular  matrix,  which  can  be  done  with  Gaussian  elimination 
and  scaling.  Next,  compute 
(26) 

=  ^2  -  hp(0  -  h’s^Oj 

J 

4-  4-  4-  4-  4- 

+  ^  ^L,0,7+l-j 

«p  [p(‘>A./c]  ^  +  gL%g|:J^o]  } 

^  =  exp  [p(0A,/c]  -f  ^(|^  +  ^(|Jr  +  cl)/]*  fo*”  >  2. 

Finally,  compute 
(27) 


’^=  /  daG[;J+(AA.-,2J/i-h2A-a)  /  dyp(0(3  -  a')  /  dV' Gj(+^’"(0,  a' -  a«)GL‘^(/‘^.-.  •»") 

'Jo  Jo  Jo 

p(0(2JA  +  2A  -  0  +  G^^\hXi,2Jh  +  2h-  3)0 

O  /  dj'G(^+0-(o,,_3/)G(;)+(AAi,j')+  /  d5p(‘)(2JA-0° 

Jo  Jo 

o  j^ds'  g(^-^O-(0,  ,  _  y)  [G(0+(;,Ai.  5'  +  2h)  +  AGW-;,).iG[:'^+(AA,,  y)]'j 
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(28) 

eg  =  f  { [l  +  ACg+i)  o]  +  2A  -  5)  +  M‘)(2/i)Gg‘J-(0. 2Jk  -  5)}  x 

X  G[;^+(Ui,  5)  +  [7j(0(2J/i  +  2/i  -  5)  +  /iGg+(.j_op(‘)(2j;i  +  2h-  j)+ 

+ hGi';;^(i).ip^<){2jh  - ,)]  Gii+^)-(o. ,)  j + 

+  [1  +  f^G<^Mi).o]  ^ 

X  r^d3  [jc(0(2j;i  +  2h-  3)p^'){3)  +  G'^+*’"(0,  2Jh  +  2k-  5)G|;^"(0,  3)] 

J2h  ^  ^  \ 

d;i  =  '‘jeSj;;-  ['i“>(2-"‘) + ciiSo..] + gKJ-  [■!'‘>(2'>) + c%..)  |- 
+ k’  (cSS"  [cSS'clii:. + c'JiT.ci;.;;)' + c!;S'Gi;.fct.] + 

1 


+ cl:, [2c!;j;;'i'‘>(2-"‘) + 2c!;«>-,(<)(2).) + 

+ gS;’;o..{°‘i!.w‘  [-j1‘1(2  a)+ ^■'1(2^'>)) + p1‘H2J'-2).)ok;;- + oKJi.ichj).)}] 


+  0{h^) 


Of  all  the  terms,  it  is  Co| j  that  uses  the  most  ink,  but,  because  it  has  no  integrals,  Co|j  is 
easily  computed.  It  is  the  triple-integral  term  in  (27)  that  is  the  most  computationally 
intensive  of  the  terms. 

It  was  already  noted,  in  the  second  paragraph  of  this  section,  that  the  long  equations 
in  this  section  contain  no  new  ideas.  The  new  ideas  in  this  paper  are  all  in  section  1  and 
in  the  paragraph  that  includes  equation  (10), 

4  Related  Work 

A  time-domain  inverse  problem  for  heterogeneous  dispersive  media  was  solved  previously 
in  the  setting  of  horizontal-shear  acoustic  waves  [12].  The  essential  difference  between  that 
inversion  method  and  the  one  presented  here  is  the  way  in  which  some  crucial  matrix 
entries  are  computed.  In  the  present  paper,  the  entries  of  the  matrix  M  and  the  vector 
[vi,  V2, . . . ,  v^r]  in  (12)  are  computed  directly,  as  described  in  (13)'-(14)  and  in  section  3.  One 
of  the  most  interesting  new  ideas  in  [12],  however,  is  the  method  by  which  its  analogous  ma¬ 
trix  entries  are  computed;  such  a  computational  method  would  later  be  called  “an  implicit 
scheme”  [l3].  The  implicit  scheme  in  [12]  uses  numerical  experiments  to  determine  matrix 
elements  in  a  sensible  way  that  is  validated  numerically  in  [12].  The  method  in  the  present 
paper,  however,  uses  direct  computation  in  place  of  numerical  experiments.  In  short,  [l2] 
has  the  first  implicit  scheme  for  heterogeneous  dispersive  (acoustic)  media,  and  the  present 
paper  has  the  first  explicit  scheme  for  heterogeneous  dispersive  (electromagnetic)  media. 
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5  Conclusion 

An  inverse  scattering  problem  for  heterogeneous  dispersive  media  is  defined  in  this  paper’s 
first  paragraph.  The  statement  of  the  problem  is  simplified  using  the  new  gauge  ideas  in 
section  1.  The  key  to  understanding  the  layer-stripping  idea  at  the  center  of  this  paper’s 
algorithmic  solution  to  the  inverse  problem  is  to  first  understand  how  the  derivation  in  (10) 
follows  from  equation  (9). 
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Propagate  ,  n  and  Related  Topics,  Including  inversion 
1.  INTRODUCTION 

The  computation  of  electromagnetic  pulses  in  dispersive  media  is  a  highly  devel¬ 
oped  field.  For  instance,  a  single  paper, ^  published  in  1976,  contains  numerics  for 
the  propagation  of  TE-  and  TM-polarized  electromagnetic  pulses  that  are  incident 
obliquely  on  an  inhomogeneous,  anomalously  dispersive  medium.  Computational 
electromagnetics  has  developed  so  extensively  since  1976  that  it  now  appears  that, 
given  enough  computer  resources,  one  can  compute  the  propagation  of  just  about 
any  single  pulse  through  just  about  any  single  medium.  But  these  studies  of  single 
pulses,  even  of  millions  of^ingle  pulses,  have  not  demonstrated  that  tvery  microwave 
pulse  travels  through  water^  with  one- ninth  the  speed  of  light  in  vacuum.  That  fun¬ 
damental  factor-of-nine  effect,  which  appears  to  have  been  unnoticed  until  now,  is 
established  here  by  studying  an  electromagnetic  wave  equation  and  its  scattering 
operators,  which  are  the  natural  places  to  find  broadly  applicable  rules,  that  govern 
propagation.  We  will  show  that  the  time  of  arrival  of  transmitted  pulses  in  anoma¬ 
lously  dispersive  media  is  related  to  a  slow  speed  given  by  the  DC  phase  velocity  in 
each  medium.  This  paper  has  several  other  new  results,  which  relate  to  the  widths 
and  peak  amplitudes  of  pulses,  and  to  quantities  that  resemble  power  density.  These 
new  results,  concerning  broad  classes  of  pulses,  are  vahdated  here  using  standard 
numerical  methods;  a  new  numerical  method  for  estimating  errors  is  also  developed 
and  used,  and  some  results  of  laboratory  experiments  on  pulse  propagation  in  a 
muscle-equivalent  matericJ  are  explained.  Our  results  will  be  shown  to  be  helpful  in 
proposing  optimum  sample  lengths  to  be  used  in  Time  Domain  Spectroscopy  studies 
for  the  accurate  determination  of  the  infinite-frequency  and  static  permittivities  of 
Debye-type  dispersive  media.^“^  Also,  as  shown  in  Section  4,  our  results  can  form 
the  bcLsis  for  a  sensitivity  cinalysis  of  the  dependence  of  the  medium  response  on 
the  parameters  obtained  from  different  fits  to  the  same  band-limited  experimental 
data. 

This  work  was  done,  in  part,  to  assist  in  the  development  of  health-and- 
safety  regulations  for  electromagnetic  pulses  in  human  environments.  Our  goal 
was  to  develop  methods  to  support  the  regulation  of  basic  quantities  such  as  the 
peak  amplitude  and  the  power  density  of  incident  pulses,  so  that  it  woiild  not 
be  necessary  to  regulate  every  detail  of  a  pulse’s  time  trace.  Toward  that  end 
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we  formulated  a  problem  in  making  inferences  from  incomplete  data.  We  asked: 
Knowing  only  the  peak  amplitude  and  power  density  of  an  incident  pulse,  what  can 
be  said  about  the  peak  amplitude  of  the  propagated  pulse?  We  aJso  asked  what 
the  incomplete  data  would  imply  about  the  values  inside  the  dispersive  medium 
of  the  time  derivative  of  the  magnetic  field  H  and  of  a  different  quantity  that  is 
related  to  power  density.  Of  those  three  quantities — peak  amplitude  and  power 
density  and  dtH — it  is  dtH  whose  size  is  most  closely  linked  to  the  time  scale  of 
short-risetime  pulses;  further,  dtH  is  particularly  important  because  it  would  be 
largely  responsible  for  elrfctromotive-force  currents  in  any  circuit-like  structure  that 
is  inside  a  dispersive  medium.  Our  results  in  this  matter  of  incomplete  data  are 
quite  concrete.  We  will  show,  for  instance,  that  the  peak  amplitude  of  the  electric- 
field  part  of  a  propagated  microwave  pulse  is  zdways  less  than  0.150  V/m,  for  depths 
greater  than  2.00  mm  in  water,  whenever  the  incident  electric  pulse’s  pe’ak  amplitude 
is  less  than  1.00  V/m  and  its  power  density  is  less  than  5.29  x  10“^^Watt/m^, 
regardless  of  the  other  details  of  the  microwave  pulse’s  time  trace.  We  have  similar 
results  for  power  density  and  for  dtH.  The  development  of  such  general  rules  for 
pulse  propagation  may  put  the  computationzJ  basis  for  pulse-safety  standards  on  as 
firm  a  basis  as  for  the  existing  standards®  for  continuous  waves  and  periodic  wave 
trains. 

This  paper’s  dispersion  models  and  time  scales  are  motivated  by  laboratory 
experiments.  We  use  a  one-term  Debye®  model  that  fits  laboratory  data  for  water^ 
up  to  100  GHz,  and  our  numerical  tests  involve  pulses  with  or  without  DC-frequency 
content  whose  time  scedes  are  characteristic  of  short-pulse  radar.  Our  methods 
apply  to  all  other  Debye-like  media,  and  can  be  generalized  for  the  two-term  and 
five-term  Debye  models  that  fit  laboratory  data  for  muscle  and  muscle-equivalent 
materials.®"^  Our  methods  can  also  be  generalized  for  non-Debye  media. 
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2.  PDE  FORMULATION  AND  RESULTS 


A.  PDEs 

The  equations  governing  the  scattering  and  propagation  of  an  obliquely  incident 
piilse  on  a  homogeneous  dispersive  hadf-space  occupying  z  >  0  are  the  time-domain 
Maxwell’s  equations  for  the  fields  Ey.  This  set  of  equations  is  coupled  through 

a  polarization  current  to  an  ordinary  differentiad  equation  that  describes  the  evo¬ 
lution  of  an  orientational  polarization  {Py)  mechanism  of  Debye  type®  (a  relaxation 
process):  r^  +  Py  =  AeEy,  where  Ae  =  e, -£oo,  e#  and  £„  are  respectively  the  zero- 
and  infinite-frequency  relative  permittivities,  and  r  is  the  dielectric  relaxation  time. 
This  o.d.e.  together  with  the  constitutive  law,  Dy  =  Coi^coBy+Py),  result  in  the  model 
frequency-domadn  relative  permittivity  e(u;)  =  goo  +  i  where  £<,  is  the  permittivity 
of  vacuum.  This  model  is  fitted  to  frequency-domain  experimental  data  for  a  range  of 
frequencies  u)  in  order  to  fix  the  Vcirious  medium  pcirameters.  Typical  values  for  water 
in  the  microwave  frequency  range  cire  e,  =  80.35,  £00  =  1-00,  r  =  8.13  psec.  The 
phaise  velocity  of  each  frequency  component  in  such  a  medium  is 
with  c  being  the  speed  of  light  in  vacuum.  In  the  subsequent  analysis  and 

v^^**(oo)  will  arise.  Finadly,  operationzd  considerations  fix  the  pulse  shape,  /,  and 
its  duration,  Tp. 

The  electric  field  incident  on  the  half-space  from  the  air  side  (z  <  0)  is  a 
plane  pulse  z,  t)  =  f{t  —  x  sin  ^.-nc/c  —  ^  cos  ^.-nc/c)  of  duration  Tp.  We  assume 

the  pulse  hcis  been  in  contact  with  the  interface  since  t  =  —00.  On  the  interface, 
z  =  0,  the  total  electric  field  is  Ey{x,Q,t)  =  g{t  —  xjv)  where  v  =  c/sin^inci  tte 
totzd  field  is  known  by  direct  measurement  of  either  the  field  on  the  interface  or  of 
the  scattered  field  in  z  <  0.  Defining  the  time-like  Vciriable  ^  t  —  x/u  we  find 
that  Ey(x,z,t)  =  E„(0,z,O  By(2,0,  ffx.x(x,z,t)  =  ff.,40,2,0  and 

that  =  —Ey(z,^).  Changing  coordinates  (x,z,t)  -)•  (z,^)  in  the  resulting 

one-dimensional  system,  and  eliminatihg  Hy  through  differentiation  with  respect  to  ^ 
and  Py  by  using  the  operator  i,  we  obtain  a  single  third-order  partial  differential 
equation  for  the  electric  field  E  —  Ey  (shown  in  factored  form), 

-  cM{d^  +  Cod,)E  +  -  ci5,)(a«  -b  cxd,)E  =  0,  z  >  0,  (2.1) 

T 

where  =  c/(y^cos  0  -  l-bA£/(eoo  cos^  <^.„c),  and  ci  =  Co/^/^.  For  <t)iru:  =  0, 
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C  =  t.  The  signaling  problem  for  (2.1)  is  completed  by  giving  the  boundciry  condition 
E(0,^)  =  g(^),  and  the  initial  data  ^(z,0)  =  E^[z,0)  =  =  0.  Hi  follows 

once  E  is  known,  and  H,.  =  -^  The  following  results  are  derived  in 

Subsection  3A. 

Equation  (2.1)  describes  the  propagation  of  all  possible  waves  of  different  or¬ 
ders,  and  their  corresponding  speeds,  that  can  be  excited  by  an  eirbitrary  pulse.  The 
coefficients  exhibit  an  explicit  dependence  on  the  angle  of  incidence,  and  on  the  pa¬ 
rameters  that  describe  the  medium.  It  is  a  strictly  hyperbolic^  partial  differential 
equation  since  the  prindipad  part  of  the  operator  has  real  distinct  eigenvalues  (three 
eigenvalues,  ±c<,  and  0),  and  a  complete  set  of  eigenvectors.  Causality  follows  from 
this  last  sentence.  The  characteristic  contributed  by  the  zero  eigenvalue  can  be  visu- 
adized  by  considering  that  d(  +  C2di  =  d(  when  C2  =  0.  The  mziin  feature  of  (2.1)  is  the 
two  wave  equations  exhibiting  distinct  speeds,  Cg  and  Ci.  Pulse  propagation  is  gov¬ 
erned  by  these  two  speeds  in  mutuedly  exclusive  spatial  regions.  Disturbances  mmnly 
described  by  the  principed  part  of  the  operator  in  (2.1)  will  be  called  high-order  waves, 
while  those  described  by  the  remzdning  operator  will  be  called  lower-order  waves.*  The 
speed  Cl,  while  not  a  chriracteristic  speed  (it  is  sub-characteristic,  Ci  <  Co),  is  impor¬ 
tant  in  the  analysis  and  has  several  ramifications  for  experiments.  Also,  ci  =  v^^**(0) 
and  Co  =  v^^'*(co);  i.e.,  the  meun  distmbances  will  propagate  with  the  distinct  speeds 
which  are  equal  to  limiting  values  of  the  phase  velocity.  It  is  worthwhile  to  empha¬ 
size  that  experimental  data  indicates  Ci  co,  e.g.,  Ci  ~  0.1116co  for  water  in  the 
microwave  range;  the  problem  is  stiff  so  the  pulse  travels  in  the  half-space  with  either 
of  two  speeds  that  are  disparate. 

The  high-order  term  describes  the  dominant  behavior  for  depths  z  <  O(cot) 
m,  and  the  effect  of  the  lower-order  term  on  the  the  high-order  waves  is  an  exponential 
decay  with  z.  The  penetrating  pulse  propagates  with  speed  Co  in  this  shcdlow  depth 
(~  10"'*  m  for  water)  which  we  name  the  “skin-depth”  for  pulses  since  it  is  remi¬ 
niscent  of  the  well  known  frequency-domain  concept.^  From  experimentally  obtziined 
data  typical  of  tissue  r  =  0(10"*^)  sec,  and  0  =  0(10).  Thus  f  is  large,  and  we 
expect  the  bulk  of  the  penetrated  pulse  to  travel  with  the  speed  ci  since  (2.1)  is  then 
approximately  =  0.  The  main  disturbance  will  be  a  lower-order  wave. 

The  effect  of  the  high-order  term  on  the  lower-order  waves,  which  travel  with  speed 
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Cl,  is  diffusive  in  character  and  important  for  z  >  0(c<,r)  m.  The  mciin  response 
diffuses  around  the  ray  z,  =  ci^  on  which  the  peak  of  the  response  is  found.  The 
peak  amplitude  on  the  sub-characteristic  ray  decays  as  or  as  l/r/f  (for  fixed 

depth).  A  consequence  of  this  is  that  the  peak  of  the  energy-like  quantity  will 
decay  as  1/z  (1/^). 

The  response  will  also  depend  crucially  on  the  pulse  duration.  This  parameter 
appears  through  scaling  ^  with  Tp,  and  z  with  cTp.  Now  ^  — )•  and  Co  and  ci  are 

normalized  by  c.  Pulses  with  appreciable  amplitude  most  often  have  Tp  ~  10~®  — 10“*° 
sec,  so  is  still  large.*  Pulses  that  zire  long  with  respect  to  the  relcixation  time 

(Tp  ->  co,  or  equivalently  if  r  — 0)  will  propagate  unattenuated  in  the  half-space 
with  amplitude  equal  to  the  DC  Vcdue  of  the  frequency-domain  transmission  coefficient 
regardless  of  the  pulse’s  DC-frequency  content.  The  field  just  after  the  interface  (no 
“skin-depth”  since  CqT  — >  0)  satisfies  a  lossless  wave  equation  with  speed  Ci.  On  the 
other  hand,  very  short  pulses  (Tp  0,  or  equivalently  if  r  — >  co)  wiU  not  penetrate 
far.  In  this  case  the  electric  field  in  the  half-space  (since  now  CqT  — >  co)  sees  a  high 
constemt  conductivity  medium  thus  it  satisfies  a  telegraphers  wave  equation  whose 
far-field  is  the  diffusion  equation. 
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B.  Green  Functions 

This  subsection  describes  some  rules  of  wave  propagation  that  are  derived  from 
time-domriin  Green  functions.  The  history  of  these  Green  functions  is  reviewed  in 
Subsection  3B.  We  will  first  state  some  results  involving  upper  bounds  on  prop¬ 
agated  peak  amplitudes  and  power-density-type  quantities.  These  upper  bounds 
are  easily  computed  and  they  are  independent  of  the  detailed  nature  of  the  incident 
fields.  The  bounds  are  developed  for  normal  incidence  here  and  for  oblique  incidence 
in  Subsection  3B.  The  present  section  concludes  with  a  description  of  wave  speeds 
and  with  brief-pulse  an<f  long-pulse  approximations.  These  rules  are  all  illustrated 
numericzilly  using  the  water  parameters  e,  =  80.35,  Coo  =  1.00,  and  r  =  8.13  psec 
from  Subsection  2A;  and  the  results  are  easily  generalized  to  other  Debye  media 
and  to  non-Debye  media.  The  necessary  derivations  and  numerical  validations  are 
in  Subsection  3B  and  Appendix  A. 

For  normal  incidence,  let  the  y-polarized  incident  electric  field  be  /(t  —  z/cq) 
in  the  air-fiUed  hcdf-space  z  <  0.  In  the  water-filled  half-space  z  >  0,  the  resulting 
electric  field  is 


E(z,t)  =  E,(z.t)  =  exp  (g  X  10-^m)  ^  -  ^). 

(2.2) 

The  Green  function  Ge('5,  t)  is  graphed  in  Fig.  1  for  several  depths  z  >  0.  For  the 
boundary  z  =  0,  Ref.  9  derives  GE(0,t)  =  —t~^  exp[— 1/(2.00  x  10”^^s)]/i[t/(2.05  x 
10~^^s)],  where  1\  is  the  modified  Bessel  function  of  the  first  kind;  the  oblique- 

incidence  generalization  is  (3.12).  The  magnetic  field  ^r3,(z,i)  =  —  zEy{z,  a)! hq 

also  has  a  Green-function  representation  similar  to  (2.2). 

Many  of  our  new  results  are  based  on  (2.2)  and  Fig.  1.  Safety  standards 
may  be  affected  by  this  type  of  analysis;  so,  in  Appendix  A,  we  show  how  to 
estimate  the  percentage  error  in  computations,  with  the  following  results  for  our 
computations:  (1)  The  pointwise  numerical  error  in  Ge(0.500  mm,  t)  is  no  more 
than  1.70%  of  the  peak  value  (with  respect  to  t)  of  |Ge(0.500  mm,  t)|;  (2)  The 
pointwise  numerical  error  in  Ge(4.00  mm,  t)  is  no  more  than  0.800%  of  the  peak 
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Fig.  1.  The  time  dependence  of  the  Green 
depths  2. 
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action -Ge  (.2,4)  for  water  at  several 
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value  of  |C?e(4.00  mm,  £)|;  and  (3)  For  intermediate  depths,  the  relative  error  in 
G^{z^t)  decreases  monotonicaUy  from  1.70%  at  z  =  0.500  mm  to  0.800%  at  z  = 
4.00  mm.  The  closed-form  expression  for  Ge(0,£)  is  exact.® 

We  will  use  the  following  three  norms: 


/  dt\h{z,t)\ 
Jo 


il^(^1-)ll2  =  (  /  cl£|h(z,£)|^ 


1/2 


i/i(z,  *)lIoo  =  least  upper  bound  of  |A(z,-)|, 


(2.3) 


where,  for  each  depth  z,  the  least  upper  bound  ||h(z,*)||oo  is  evaluated  with  respect 
to  i.  Then,  for  each  depth  z  from  0.500  mm  through  4.00  mm. 


<  l/(-)looexp 


(0.202)  ||/(•)||oo,] 

F2{z)\\f{-)h  J 


(2.4) 


regardless  of  the  detailed  nature  of  the  incident  electric  field  /(£).  The  right  side 
of  (2.4)  is  easily  computed,  given  the  functions  Fi{z)  =  |Ge(2:, •)|2  -^2(2:)  = 

j(?E(2;,  •)!<»>  which  are  graphed  in  Fig.  2.  Inequality  (2.4)  defines  upper  bounds  on 
the  peak  amplitudes  of  JS(z,  •).  This  inequality,  and  all  of  our  other  Green-function 
results,  are  validated  numerically  in  Subsection  3B.  The  upper  bound  on  the  right 
side  of  (2.4)  is  almost  attained  in  one  of  the  numerical  vzdidations.  In  that  sense, 
the  upper  bound  is  sharp. 

We  wiU  now  show  how  relation  (2.4)  could  be  used  in  a  safety  standard  for 
the  peak  amplitudes  of  internal  electric  fields.  Suppose,  for  this  hypothetical  exam¬ 
ple,  that  it  has  been  determined  that  peak  electric  fields  must  be  no  greater  than 
0.200  V/m  at  depths  greater  than  1.00  mm,  and  no  greater  than  0.150  V/m  at 
depths  greater  than  3.00  mm.  That  hypothetical  internal-field  standard  is  trans¬ 
lated,  using  (2.4)  and  Fig.  2,  into  a  more  easily  regulated  standard  on  incident 
electric  fields  /(£).  The  easily  regulated  standard  is 
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(X  10"^  sec 
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0.0005 


0.0015 


0.0025 


0.0035 


Depth  (meter) 
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1 

0.5 
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Rg.  2.  Th.  L,  norm  F,{z)  nnd  the  L„  norm  ^i(x)  of  Iho  Green  funetion  Ge{z,  ■)  for 

water.. These  norms  and  Eqn.  (2.4)  reduce  the  upper-bound  computations 
to  a  calculator  exercise. 
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(X  10'Vsec) 


(1)  Peak  value  of  |/(i)|  <  0.740  V/m  or 


(2) 

(3) 


Peak  value  of  |/(f)|  <  40,  000  V/m  and 
Peak  value  of  \f{t)\  <  50,  000  V/m  and 


or 


(2.5) 

By  reading  the  two  graphs  in  Fig.  2  and  using  a  calculator,  one  can  use  (2.4)  to 
show  that  any  incident  field  that  satisfies  item  (1)  or  item  (2)  or  item  (3)  of  (2.5) 
is  guaranteed  to  produce  internal  fields  that  comply  with  the  hypothetical  internal- 
field  standard,  regardless  of  all  other  details  of  f{t). 

Upper  bounds  also  exist  for  quantities  that  resemble  power  densities.^®  In 
particular,  for  any  incident  field  /(/)  and  for  each  depth  z  from  0.500  mm  through 

4.00  mm,  the  power-density-type  term^°  [||JS(z,  •)||2]^  =  f^dt\E{z,t)\^  satisfies 


i  ^  “Ke-isxio-O 


(0.202)  |/{•)b,11“ 


(2.6) 

Subsection  3B  numerically  vadidates  the  inequality  in  (2.6),  showing  that  the  upper 
bound  on  the  right  side  of  (2.6)  is  almost  attained  in  at  least  one  case. 

We  win  now  show  how  the  upper  bounds  in  (2.6)  could  be  used  in  a  safety 
standard  for  an  power-density-type  quantity  related  to  internal  electric  fields.  Sup¬ 
pose,  for  this  hypothetical  example,  that  it  has  been  determined  that  the  power- 
density-type  quantity  J*Q*”dt|Fl(z,  t)p  must  be  no  greater  than  1.50  x  10~^^  V^s/m^ 
at  depths  greater  than  1.00  mm,  and  no  greater  than  1.00  x  10“^^  V^s/m^  at  depths 
greater  than  3.00  mm.  Equation  (2.6)  and  Fig.  2  translate  this  hypothetical  stan¬ 
dard  for  internal  fields  into  a  more  easily  regulated  standard  on  incident  electric 
fields  f{t): 
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m  or 


(1)  /  dt\f{t)\^  <  2.40  X  lO-^V^s/ 

Jo 

(2)  dt\f[t)\^  <  36.0  V^s/m^  and  J  di|/(f)|  <  3.90  x  10“^^  Vs/m  . 


(2.7) 


Subsection  3B  shows  that  some  pulses  almost  attain  the  upper  bounds  in  (2.6), 

which  was  used  to  obtain  (2.7). 

« 

The  upper-bound  concepts  in  (2.4)  and  (2.6)  are  easily  extended  to  mag¬ 
netic  fields  and  their  time  derivatives,  and  to  oblique  incidence.  Subsection  3B  has 
numerical  residts  for  all  of  those  extensions.  One  can  see  there  that  making  the 
angle  of  incidence  more  oblique  will  decrease  the  penetration  into  the  medium  of 
power-density-type  quantities  and  also  peak  electric  and  magnetic  fields.  A  related 
closed-form,  modified-Bessel-function  expression  for  the  oblique-incidence  reflection 
kernel  R^{t)  =  GE,0(O,t)  is  given  in  (3.12). 

We  will  now  state  some  Green-function  results  concerning  wave  speeds. 
These  results  axe  derived  in  Subsection  3B.  For  simplicity,  the  results  are  stated 
for  normal  incidence.  Our  first  conclusion  is  that  the  main  bulk  of  an  electromag¬ 
netic  pulse  travels  through  water  with  speed  cq  for  0.3  mm,  and  then  slows  until,  for 
all  depths  beyond  0.7  mm,  the  pulse  travels  with  the  constant  speed  co/9.0.  That 
behavior  contrasts  with  the  wavefront  speed,  which  is  mathematically  well  defined 
but  is  not  tdways  observable  in  a  laboratory.  The  wavefront  speed  is  precisely  cq 
for  all  depths. Section  4  discusses  various  Debye  models  that  are  consistent,  to 
within  about  10%,  with  the  band-limited  water  data^  used  here.  The  large-depth 
speeds  (all  «  co/9.0)  for  those  Debye  models  vary  by  only  about  10%.  The  shedlow- 
depth  speed  and  the  wavefront  speed  of  any  Debye  model,  however,  are  both  equal 
to  Co  =  If  y/floiZo-  The  shallow-depth  and  wavefront,  speeds,  consequently,  change 
considerably  as  one  varies  the  Debye  parameter  Coo  from  1.00  through  10.0,  as 
described  in  Section  4.  Therefore,  a  measurement  of  the  wavefront  speed  or  the 
shallow-depth  speed  would  determine  the  Debye  parameter  Coo* 

We  conclude  with  some  Green-function  results  concerning  short-pulse  and 
long-pulse  approximations.  The  propagation  of  any  finite-valued  incident  electric 
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pulse  f{t)  is  given  simply  by  (2.2)  and  Fig.  1.  We  get  additional  insight  by  con¬ 
sidering  approximations  for  what  we  will  call  elemental  pulses  /*:  An  element  /e(t) 
is  zero  except  on  a  single  time  interval,  during  which  it  is  either  strictly  negative 
or  strictly  positive.  For  instance,  a  square  pulse  is  an  element,  but  a  one-cycle 
sinusoid  is  not  an  element.  Elements  are  important  because  any  incident  pulse  / 
is  a  sum  of  positive-valued  elements  and  negative-valued  elements.  If  the  dura¬ 
tion  of  an  element  f^{t)  is  much  briefer  than  30  psec,  then  the  propagated  pulse 
element  is  approximately  ll/e(‘)lli^E(-2^>0  (see  Fig*  1)  depths  greater  than 

0.7  mm.  If  the  duration  of  an  element  /e(t)  is  much  longer  than  50  psec,  then  the 
propagated  pulse  element  is  approximately  0.2/e  [t  —  9.0(z  —  1  mm)/c  -f  17  psec]  = 
0.2/e  {t  -  9.0z/c  -  13  psec)  for  aU  depths  from  1  mm  through  the  depth  at  which 
the  duration  of  GE{z,t)  becomes  comparable  to  the  duration  of  /e(0-  These  ap¬ 
proximations  are  derived  in  Subsection  3B. 
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3.  DERIVATIONS 


A.  PDEs 

To  extract  from  (2.1)  the  equation  describing  the  ezirly-time  evolution  (in  the  “skin- 
layer”)  of  the  response  we  set  everywhere  in  (2.1)  ~  —Co^z  except  in  the  operator 

di  -f  Codz  since  it  expresses  the  propagation  of  the  high-order  waves.  Any  other  terms 
in  the  resulting  equation  will  describe  the  effect  of  the  lower-order  waves.  The  main 
disturbance  for  early  times  is  modeled  by 


(5f  -|-  Codz)E  -b 


(3  cl- cl 
T  2cl 


E  =  0]  Z  <  CoT, 


(3.1) 


subject  to  the  boundary  condition  5(0,^)  =  g{()’  The  solution  of  (3.1)  is 


=  aU - )e-'^P 

Co 


_  f £»/£co__l\ 

\2C0S^  ^inc)  \CoTJ 


(3.2) 


We  see  that  the  response  decays  exponenticilly  in  a  thin  region  of  depth  z  ~  O(cor), 
where  the  speed  of  propagation  is  Cg.  Note  that  the  decay  constant  is  inversely 
proportioned  to  cos^  (f>inc  thus  normed  incidence  will  result  in  the  greatest  amplitude 
in  the  medium.  To  describe  the  evolution  of  the  lower-order  waves,  which  travel  with 
speed  ci,  we  set  in  (2.1)  ~  —cidz  except  in  the  operator  d(  -f  cidz  which  expresses 

the  hyperbolic  nature  of  the  lower-order  waves.  The  main  disturbance  is  now  modeled 

by  2  2 

{d^  +  c,dz)E  =  ^^^dlEi  z>CgT.  .  (3.3) 

The  boundary  condition  is  approximately  E(zo,  ^*)  =  ^(^*)>  where  Zg  is  the  depth  aifter 
which  (3.1)  no  longer  applies,  is  the  time  with  origin  at  Zgjcg  (the  time  it  takes 
for  the  pulse  to  reach  Zg  in  the  “skin-depth”),  and  h[^)  represents  (3.2)  evaluated 
at  Zg.  Equation  (3.3)  is  an  advection-diffusion  equation,  and  describes  the  response 
cifter  a  depth  of  0[cgT)  m.  The  peak  of  the  response  is  on  the  sub-characteristic 
ray  z,  =  ci^*.  The  solution  is  very  easily  obtziined  from  the  solution  of  the  diffusion 
equation.  It  is^^ 


X  /'  dn- 

Jo 


0  (r-«)t 


exp 


-  Ci(^'  -  ^)P 


(3.4) 
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Various  techniques  can  be  used  to  estimate  the  integral  in  (3.4)  since  (3 fr  is  large. 
Here  we  are  interested  only  in  the  primary  behavior  of  £?  as  a  function  of  depth.  For 
^  1  the  response  is 


E‘(z,i')  =  h{0h 


1 

/?  1 

[(-5  - 

\ 

1 

2r(c2  -  c?) 

1 

1 

(3.5) 


On  z,  =  we  find  that  7nax{E^}  ~  or  ~  yJJ'.  This  is  verified  with  the 

numerical  experiments  in  Section  4. 
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B.  Green  Functions 

This  subsection  derives  our  Green-function  results  in  the  order  in  which  they  are  de¬ 
scribed  in  Subsection  2B.  Equation  (2.2),  for  instance,  is  a  Green-function  represen¬ 
tation.  These  Green  functions  have  become  a  standard  technique  in  computational 
electromagnetics.  They  were  first  developed  for  non-dispersive  media,^^  and  were 
then  used  to  compute  fields  in  dispersive  media.^'^  That  dispersive- medium  work 
has  not  yet  been  published,  owing  to  the  death  of  R.  Krueger,  but  generalizations 
are  avziilable.^®"^^  The  Green  function  programs  used  in  the  present  paper  were 
developed  by  the  authors  of  Ref.  16,  along  the  lines  of  Appendix  A  of  that  paper. 
The  Loo  and  L2  norms  in  (2.3)  have  special  physical  significance.  The  L2 

norm  is  important  because  (c£o/2)  (HFIj)^  is  the  power  density,  whose  mks  units 
are  Watt/m^,  of  an  electric  pulse  in  free  space.^°  Consequently,  we  wiU.  focus  on  the 
peak-Vtilue  (p  =  00)  and  power-density  (p  =  2)  cases  of  the  inequalities 


|f:(r,.)t<e-“|/(.)|p+|G(r,.)l.l/(-)l,.  (3.6) 

which  are  obtained  by  applying  the  Young  theorem^®  to  (2.2),  and  for  which  a  = 
1.63  X  10^  m“^.  In  particular,  although  (3.6)  is  valid  whenever  1  <  p,  q,  r  <  00 
satisfy  r“^  =  p~^  +  q~^  —  1,  we  are  most  interested  in  the  cases 


ll■S(^,•)loo  <  e  “^ll/(•)ioo  +  min 


|Ge(^,-)Ii1I/(-)IIoo, 

l|GE(^,-)bi/(-)ll2, 

lGE(z,-)llco||/(-)ix 


(3.7) 


and 


(il-S(z,.)l|2]^  < 


e 


—  az 


ll/(•)ll2  +  mini 


I1Ge(^,-)IIi|/(-)1|2, 

iGB(2:,-)B2i/(-)lli 


(3.8) 


In  numerical  computations  for  water,  ||C?b(2!,  •)li  was  observed  to  decrease  slowly 
and  monotoniccdly  from  0.2019  at  0.5  mm  to  0.2014  at  3.5  mm.  It  is  as  if,  for  those 
depths  in  water,  the  advection-diffusion  equation  (3.3)  were  approximated  by  a  heat 
equation  and  GE[z,t),  which  is  positive  valued  for  depths  beyond  0.5  mm,  were 
analogous  to  a  temperature  distribution  whose  total  conserved  heat  is  proportional 


IV  -  15 


to  |Gb('2,  oil-  Equations  (3.6)-(3.8),  the  almost-constant  nature  of  ||<je(^j ')ii  >  and 
numerical  Green-function  computations  produced  the  results  in  (2.4)-(2.7). 

We  now  consider  five  numerical  examples  that  vadidate  the  inequalities  in 
(2.4).  We  will  see  that  the  minimum  upper  bounds  in  (2.4)  are  almost  attained  for 
some  incident  pulses.  For  these  five  examples,  we  choose  the  following  hypothetical 
restrictions  on  the  incident  electric  pulses  f{t): 

l/(Oico<1.00  V/m  and 

# 

l/(-)ll2  <  6.32  X  10"®  Vs^/Vm  and  (3.9) 

|/(•)il  <  4.00  X  10“^^  Vs/m. 

Then  (2.4)  shows  that  any  incident  pulse  /  that  satisfies  (3.9)  wiU  produce  an 
internal  field  whose  peak  amplitude  satisfies 


0.202  V/m, 

fi  V  in-5m  ) 

-f  min 

(6.32x10-®  Ys^^^/m)Fi{z), 

\d.15  X  10  “m / 

(4.00  X  10-^^  Vs/m)E2(z) 

(3.10) 

where  Fi{z)  and  F2{z)  are  graphed  in  Fig.  2.  Each  sum  of  the  exponential  in 
(3.10)  and  a  term  from  the  “min”  clause  of  (3.10)  yields  one  of  the  three  top-most, 
boldface  curves  in  Fig.  3.  Relation  (3.10)  guarantees  that  the  depth-dependent  peak 
amplitudes  of  E{z,  •)  are  less  than  the  minimum  of  the  three  boldface  upper-bound 
graphs.  That  prediction  was  tested  using  five  incident  pulses  /  that  comply  with 
(3.9).  Those  incident  pulses  are:  (1)  a  40-psec-duration  square  pulse  with  1-V/m 
amplitude;  (2)- the  absolute  value  /2(t)  =  1/4(01  4-cycle,  80-GHz  sinusoid  in 

item  4  below;  (3)  the  absolute  value  /^{t)  =  1/5(01  1-cycle,  80-GHz  sinusoid 

in  item  5  below;  (4)  a  4-cycle  80-GHz  sinusoid  with  1-V/m  amplitude;  and  (5)  a 
1-cycle  80-GHz  sinusoid  with  1-V/m  amplitude.  The  norms  of  the  pulses,  which 
are  tabulated  below,  all  satisfy  (3.9). 
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Peak  Value  (Volt/meter) 


Fig.  3.  Five  numerical  validations  of  the  upper-bound  concept  for  the  depth-dependent 
peak  amplitudes  of  E{2,  •).  The  boldface  curves  are  upper  bounds  from  re¬ 
lation  (3.10). 


Table  1.  Five  incident  pulses  that  satisfy  the  conditions  in  (3.9).  The  first  incident 
pulse  is  1  V/m  for  0  <  f  <  40  psec,  and  it  is  0  for  all  other  times.  The  second  pulse 
is  the  absolute  vzJue  of  (IV/m)  sin[27rt/(8  X  10^°s)]  for  0  <  t  <  50  psec,  and  it  is  0 
for  all  other  times. 


Example 

Duration 

Type 

1 

40-psec  * 

square  pulse 

2 

4  cycles 

|80-GHz  sine 

3 

1  cycle 

l80-GHz  sine 

4 

4  cycles 

80-GHz  sine 

5 

1  cycle 

80-GHz  sine 

1-IU 

1  •  lb 

11 -ii 

(V/m) 

(Vs^/^/m) 

(Vs/m) 

1.00 

6.32x10“® 

4.00x10““ 

1.00 

5.00x10“® 

3.18x10““ 

1.00 

2.50x10“® 

7.96x10““ 

1.00 

5.00x10“® 

3.18x10““ 

1.00 

2.50x10“® 

7.96x10““ 

The  peak  amplitudes  of  the  five  internal  fields  JS(z,*),  corresponding  to  the  above- 
tabulated  incident  fields,  are  also  graphed  in  Fig.  3;  the  curves  for  Examples  4 
and  5  almost  overlap.  Those  peak  amplitudes  are  all  less  than  the  (boldface)  upper 
bounds  described  earlier.  The  five  examples,  therefore,  numerically  validate  the 
upper  bound  concept  in  (2.4).  Fig.  3  also  shows  that  the  upper  bounds  are  sharp 
in  the  sense  that  the  peak  amplitude  of  one  pulse  (Example  1)  aJmost  attaiins 
the  minimum  upper  bound.  That  example  involves  a  pulse  with  a  nonzero  DC- 
component^®"^^  dtf(t).  It  makes  intuitive  sense  that  the  presence  of  a  DC 
component  in  a  pulse  would  tend  to  diminish  the  attenuation  of  the  pulse  in  any 
medium,  as  an  elementary  anadysis^^  affirms  for  a  non-Debye  medium. 

Having  just  validated  the  upper-bound  concept  (2.4)  for  peaks,  we  now  val¬ 
idate  (2.6):  The  two  top-most,  boldface  curves  in  Fig.  4  correspond  to  the  upper 

bounds  (e““''|l/i2  +  |GE||ii/ll2)^  and  (e"“*ll/||2 +  iGEi2i/l|i)^  in  (2.6),  subject  to  the 
hypothetical  restriction  (3.9).  The  other  five  curves  represent  the  power-density- 

type  quantities  [||E3;(2:,  •)||2]^  produced  by  to  the  five  incident  fields  in  Table  1.  We 
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see  that  the  field  produced  by  Example  1  of  Table  1  almost  attains  the  minimum 
upper  bound  in  Fig.  4.  This  completes  our  validation  of  the  upper-bound  concepts. 

As  explained  above  (2.1),  oblique  incidence  is  taken  into  account  using  a 
simple,  widely-known  transformation  of  variables.  Using  the  transformation,  we 
obtain  numerical  results  for  dS^-incident  electric  fields  f[t  —  (i  -b  z)/(v^c)].  In 
water  (z  >  0)  the  y  component  of  the  electric  field  is 

t-z/{V2c) 

ds/(s)GE  ,450  (z,  t  “  3). 

(3.11) 

The  function  Ge,45*  {z,  t)  is  graphed  in  Fig.  5  for  several  depths  z  >  0.  The  numeri¬ 
cal  errors  in  Fig.  5  were  quantitatively  estimated  using  the  method  of  Appendix  A. 
The  results  are:  (1)  For  each  depth  z,  from  0.160  mm  through  2.88  mm,  the  er¬ 
ror  in  the  computed  values  of  GE,45»(z,t)  are  no  more  than  3.19%  of  the  peak 
value,  with  respect  to  time,  of  the  actucd  values  |GB{2)t)|;  and  (2)  The  relative 
errors  decrease,  but  not  necessarily  monotonically,  from  3.19%  at  0.160  mm  to 
2.07%  at  2.88  mm.  At  the  boundary  z  =  0  and  for  all  t  >  0,  GB,4s«(0,t)  = 
— exp[— i/(1.01  X  10“^^s)]Ji[t/(1.03  x  10“^^s)],  where  Ii  is  the  modified  Bessel 
function  of  the  first  kind.  More  generally,  the  oblique-incidence  transformation  and 
Ref.  9  imply,  for  aU  t  >  0,  that 


=  "^P(4.3ox~i0-=m) 


—  GE,0(O,t)  =  exp 


2eoocos2  0)  ^  ^'(2eoocos2  5)  •  • 


That  exact  result  uses  the  modified  Bessel  function  of  the  first  kind  to  represent 
the  reflections  in  the  cur-filled  hadf-space  (z  <  0)  that  are  caused  by  waves  that  are 
incident  obliquely  on  the  Debye  half-space  (z  >  0)  defined  by 


D( r  i\  —  j 

’  “  \£co£o-E(2,0 


z  <  0 
z  >  0  ■ 
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(3.13) 


Time  Integral  (Volt^sec/m^) 


Fig.  4.  Five  numerical  validations  of  the  upprer-bouird  concept  for  the  power- density- 
type  time  integral- /“df|F7(z,t)|2.  The  boldface  upper  bounds  are  from  rela- 


G6.45"  (X  10’Vsec) 


Fig.  5.  The  time  dependence  of  the  Green  function  Gb,45»  (^,  t)  for  several  depths  z. 
This  Green  function  is  used  for  pulses  that  axe  incident  on  water  at  an  angle 
of  45*. 
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In  particular,  the  reflected  field  is  exactly  ds/{s)iE^(t  —  j)  for  any  incident  field 
/.  Equation  (3.12)  was  used  numerically  to  vcdidate  the  z  =  0  boundary  values  of 
the  oblique-incidence  Green-function  computations  that  yielded  Figs.  5,  6  and  8. 
Equation  (3.12)  was  also  used  to  validate,  at  z  =  0,  the  previous  percentage-error 
estimates  for  oblique  incidence.  At  z  =  0,  the  estimated  relative  error  was  1.91%; 
the  true  relative  error  was  1.87%. 

Applying  the  Young-theorem  result  (3.6)  to  the  oblique-incidence  represen¬ 
tation  (3.11)  yields  obvjous  oblique-incidence  generalizations  of  the  upper-bound 
results  (2.4)-(2.6);  for  instance,  |E(0,z,-)||2  <  1/(*)|2  exp[— z/(4.35  x  10“®m)]-|- 
inin[|GB,45‘(2,-)lill/(-)ll2,  11Ge,45»(2,')I1211/(-)IIi1-  ThenormsF3(z)  =  |lGE,45»(2:,0i2 
and  ^4(2')  =  ||Ge,45*{2:,-)IIoo  are  graphed  in  Fig.  6,  and  ||Ge,45»  (2,  •)||i  was  observed 
to  decrease  slowly  and  monotonically  from  0.149  at  0.240  mm  to  0.144  at  3.00  mm. 
We  numerically  tested  these  oblique-incidence  inequeilities  using  the  five  pulses  in 
Table  1.  The  inequalities  were  validated  in  each  case,  and  the  minimum  upper 
bound  was  almost  attained  in  the  case  of  a  DC-component  pulse  (Example  1). 

We  will  now  substantiate  the  results  in  the  last  two  paragraphs  of  Subsec¬ 
tion  2B.  The  results  rely  mainly  on  (2.2)  and  Fig.  1.  Note  that  the  first  term  on 
the  right  side  of  (2.2)  represents  a  wave  that  travels  with  speed  Co  and  decays  ex- 
ponenticdly  by  a  factor  of  132  in  each  0.300  mm  interval.  Therefore  the  convolution 
term  in  (2.2)  predominates  for  depths  greater  than  0.300  mm.  The  major  features, 
such  as  the  peaks,  of  the  convolution  kernel  Ge  are  seen  in  Fig.  1  to  travel  more 
slowly  than  Cq  for  depths  greater  than  0.500  mm.  The  peak  of  GE{2,t),  for  in¬ 
stance,  is  shown  in  Fig.  7  to  travel  with  speed  cq  for  the  first  0.300  mm,  and  then 
to  slow  gradually  to  co/9.0.  (The  smciU  non-monotonic  feature  at  shcdlow  depths 
is  a  numerical  artifact  caused  by  applying  the  mcLx(-)  function  to  a  peak  that  is 
broad.)  The  numericcdly  determined  fast  speed  cq  and  the  slow  speed  co/9.0  agree 
quantitatively  with  analytical  results  in  Subsection  2A,  and  also  substantiate  the 
results  in  the  next-to-last  paragraph  of  Subsection  2B. 

The  last  paragraph  of  Subsection  2B  concerns  elemental  pulses  /e(t))  which 
are  zero  except  on  a  single  time  interval,  during  which  they  are  either  strictly  posi¬ 
tive  or  strictly  negative.  For  example,  the  Green  function  ^5(2,^)  is  an  elemental 
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^  u.uu  1  U.UUk:  0 

Depth  (meter) 

Fig.  6.  The  Li  norm  (z)  and  the  Loo  norm  Fj  (z)  of  the  Green  function  Ge,45*  (z,  •)• 

The  paragraph  below  (3.13)  relates  these  norms  to  upper  bounds  for  obliquely 
incident  pulses. 
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Time  (psec) 


Depth  (meter) 

Fig.  7.  The  time  of  arrival  of  the  peak  of  the  Green  function  at  various 

depths  in  water. 
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(X  1 000  sec 


Depth  (meter) 


Fig.  8.  The  Lj  norm  F5{z)  and  the  norm  Fi{z)  of  the  Green  function  C?h,4s*  (2,  •) 

for  magnetic  fields  produced  by  pulses  that  are  incident  on  water  at  an  angle 

of  45°.  The  last  paragraph  of  Section  3  relates  this  figure  to  upper  bounds 
for  dtH. 
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pulse  for  all  of  the  depths  graphed  in  Fig.  1.  The  statements  in  Subsection  2B  that 
concern  the  propagation  of  elemental  pulses  all  rely  on  the  following  approxima¬ 
tion:  If  an  elemental  pulse  /brief(0  is  much  briefer  than  an  elemented  pulse 

and  if  /brief(0  is  concentrated  at  the  time  foj  then  l/^dj/briefCt  ~  •*)/long("»)|  = 

l/odj/iong(t  -  •s)/bnef(-s)|  «  i/brief (Oil l/iong(f  -  k)].  That  approximation  is  asso¬ 
ciated  with  the  concept  of  (J-sequence  functions  in  the  elementary  theory  of  Dirac 
delta  functions.  We  also  note  that  the  exponential  term  in  (2.2)  quickly  becomes 
negligible  because  it  detays  by  a  factor  of  132  across  each  0.300  mm  interval  of 
depth.  These  three  results — the  negligible  exponential,  the  approximation  of  the 
convolution,  and  the  observation  that  Ge  is  an  elemental  pulse — together  with 
(2.2)  and  Fig.  1  yield  the  brief-pulse  result  in  Subsection  2B,  which  has  the  term 
||/e||iGE(2,t).  The  long-pulse  result  in  that  subsection  uses  the  additionzd  obser¬ 
vation,  below  (3.8),  that  ||Ge(2>*)IIi  approximately  constant,  and  the  time  shift 
in  the  long-pulse  result  in  Subsection  2B  also  relies  on  Fig.  7  and  results  from  the 
previous  paragraph. 

In  a  fined  matter  we  would  like  to  make  it  clear  that  we  do  not  know  what 
are  the  medical  effects  of  isolated  pulses.  We  have,  however,  developed  meth¬ 
ods  that  are  flexible  enough  that  they  may  be  useful  once  the  medical  effects  are 
known.  Although  we  originally  developed  the  upper-bound  method  along  the  lines 
of  peak  amplitudes  and  power  densities,  following  Ref.  5,  the  method  can  eas¬ 
ily  be  extended  to,  say,  the  time  derivative  of  the  magnetic  field  Hx{z,t).  To 
demonstrate  the  extension,  we  computed  the  Green  function  Gn  in  Hx{0,z,t)  = 

—  ^e~“'*/(t  —  z/coo)  +  J^/c  d3/(t  —  3)(jH('2,'S)j  /(/^oCoo)  using  the  methods  of  Refs.  13- 

17,  and  we  validated  the  computation  as  described  in  Section  4.  The  Green-function 
representation  for  ffx  yields  the  following  analog  of  (3.6): 

M0CooldtHx(zr)lp  <  e— ||aj(.)ip  +  |/(0)||IGh(z,-)1Ip  +  I1C?h(z, •)IIr||at/(*)il„  where 

r~^  =  p~^  -t-  q~^  —  1  and  1  <  p,  g,  r  <  oo.  The  norms  ^5(2)  =  ||Gh,45* (2) OIU 
and  Fe(z)  =  |Gh,45« (2, •)l|oo  for  that  inequality  are  graphed  in  Fig.  8  for  the  45°- 
incidence,  magnetic  Green  function  Gh,45'>* 
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4.  NUMERICAL  VALIDATION  AND  LABORA¬ 
TORY  EXPERIMENTS 


We  vcilidated  the  numerical  results  obtained  with  the  Green’s  function  approach  by 
compciring  the  computed  electric  and  magnetic  fields  in  a  Debye  half-space  (Section  2, 
<f>inc  =  0)  to  those  computed  with  a  finite  difference  method. Figure  9  shows 
electric  field  time-traces  computed  at  three  spatial  locations  in  the  half-space  due  to 
a  Tp  =  40-psec  square  pulse  of  initial  amplitude  1  V/m.  We  note  excellent  agreement 
to  within  a  width  of  the  line  over  an  amplitude  scale  of  10  orders  of  magnitude.  This 
indicates  a  better  than  l-part-per-billion  agreement.  (The  smadl  difference  increases  to 
two-pen-strokes’  width  at  the  shcillowest  depth,  but  only  after  the  field  itself  decreases 
by  a  factor  one  thousand.)  Also,  the  speeds  of  the  first  cirrival  and  of  the  peak  of  the 
response  can  be  deduced  from  this  graph.  We  see  that  the  first  cirrival  occurs  with 
speed  Co  =  c,  while  the  peak  of  the  response  arrives  with  the  speed  Ci  =  0.1116c  as 
predicted  in  Section  2.  Figure  10  shows  a  comparison  of  magnetic  fields,  computed 
with  the  two  methods  at  three  depths  in  the  half-space,  for  a  square-modulated 
sinusoidal  pulse  of  Tp  —  50-psec  duration  and  Ccirrier  frequency  80  GHz.  Again 
similarly  excellent  agreement  is  noted. 

Next,  we  confirm  the  analytical  results  presented  in  Subsection  2 A  by  solving, 
with  the  Green’s  function  method,  for  the  impulse  response,  f{t  —  xfc)  =  5[t  —  x/c), 
of  the  Debye  half-space  described  there  for  (f>inc  =  0  (^  =  0-  determine  the  two 
speeds  predicted  by  (2.1)  the  peak  of  the  impulse  response  was  obtziined  from  the 
simulation  results.  The  temporal  versus  the  spatial  location  of  the  peak’s  occurrence 
is  graphed  in  Figure  7.  The  slope  of  the  paph  is  the  reciprocal  of  the  speed  of  the  peak 
of  the  response.  The  relative  unimportance  of  the  characteristic  z  =  Cgt  (equivalently, 
the  wavefront  speed)  with  respect  to  the  sub-characteristic  z  =  cit  is  immediately 
evident.  The  vzdue  of  the  slope  is  given  in  Figure  7  for  two  ranges  of  depth.  We  see 
that  for  depths  of  0{cot)  the  response  travels  with  the  speed  Cq,  and  then  slows  down 
in  another  0{cot)  interveil.  This  additional  interval  will  be  much  smcJler  or  eiltogether 
absent  for  band-limited  pulses.  Eventually,  the  pulse  travels  with  the  speed  ci  for 
the  remainder  of  the  hailf-space.  The  diffusive  chziracter  of  the  meiin  disturbance  in 
the  half-space  is  exemplified  in  Figure  11  where  the  numerically  obtained  max{E^}  is 
compcired  with  the  analyticzdly  predicted  behavior  (l/->/z)  as  a  function  of  depth.  The 
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Electric  Field  (V/m) 


Time  (psec) 


Fig.  9.  The  electric  fields,  at  three  depths  in  water,  that  result  from  a  normally 
incident  40-psec  square  pulse  of  1-V/m  amplitude.  At  the  three  depths 

graphed,  the  fields  are  precisely  0  until  5.46  psec,  8.20  psec,  and  9.60  psec 
respectively. 
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Fig.  10.  The  magnetic  fields  that  result  from  a  normally  incident  4-cycle  1-V/m- 
amplitude  square-modulated  sinusoid  of  50-psec  duration. 
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Fig.  11.  The  predicted  and  computed  decays  of  the  peak  field  as  a  function  of  depth. 
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constant  of  proportioncility  for  the  l/\/z  prediction  was  fit  to  the  data  at  the  center- 
point  (4  mm)  of  the  graph.  We  note  that  although  the  slow  speed  is  achieved  very 
early  on  (z  ~0.1  mm),  it  takes  some  time  for  the  pulse  to  stzirt  diffusing  (z  ~  2.5  mm). 
This  delay  depends  on  the  frequency  content  of  the  incident  signal,  and  will  be  shorter 
for  band-limited  pulses.  From  the  discussion  and  figures  it  is  evident  that  the  envelope 
and  duration  of  the  incident  pulse  control  the  magnitude  and  nature  of  the  response 
of  the  medium.  Any  high-frequency  carrier  component  decays  exponentially  with 
depth.  Similar  behavior  was  observed  in  numericcJ  simulations  with  square  pulses  of 
various  durations.  All  long  duration  pulses  were  found  to  travel  unattenuated  with 
amplitude  T(a;  =  0)  x  max {£?"“},  and  square  modulated  sinusoidal  pulses  of  veirious 
carrier  frequencies  and  durations  (1  to  10  cycles)  behaved  like  the  Green’s  function 
with  the  cairrier  component  exponenticJly  smziU. 

We  also  checked  the  sensitivity  of  o\ir  calculations  to  the  way  in  which  we 
modeled  the  water  data  of  Ref.  2.  Reference  25  points  out  that  many  other  data  sets 
are  available  for  water,  so  we  were  satisfied  with  any  model  that  fit  the  data  in  Ref.  2 
to  within  10%.  In  particular,  we  examined  several  Debye-model  fits  to  the  data  in 
Ref.  2  for  frequencies  that  are  below  100  GHz  and  for  which  the  data  are  also  said, 
in  the  reference,  to  be  reliable.  Our  experience  in  fitting  those  data  is  that  Coo  can 
be  taken  to  be  any  number  from  1.00  through  10.0,  and  then  values  for  the  two  other 
Debye  parameters,  e,  and  r,  can  be  found  that  fit  the  data  to  within  10%.  In  that 
sense,  the  Vcdue  of  e^o  is  somewhat  arbitrciry;  for  most  of  our  simulations  we  chose 
Eoo  =  1.00,  which  is  consistent  with  an  assumption  in  Ref.  2.  The  corresponding 
Debye  parameters  zue  Coo  =  l-06  and  —  80.35  and  r  =  8.13  x  10~^^s~^  in  the 
notation  of  Subsection  2A,  or,  equivadently,  a  =  9.76  x  lO^^s”^  and  h  =  1.23  x  lO^^s"^ 
in  the  notation  of  (3.13).  This  fit  is  referred  to  as  Model  1  in  Fig.  12;  it  is  the  fit  that 
is  used  in  most  of  the  numericzJ  computations  in  this  paper.  Another  model  that  fits 
the  water  data  to  within  10%  accuracy  is  defined  by  =  5.50  and  e,  =  78.20  and 
r.=  8.1  X  10“^^s~^,  or,  alternatively,  a  =  8.98  x  lO'^s"^  and  h  =  1.23  x  10“s"^  This 
second  model  is  used  only  in  Fig.  12,  where  it  is  called  Model  2.  The  propagation  of  an 
incident  5-cycle  8-GHz  1-V/m-amplitude  square-modulated  sinusoid  was  computed 
for  these  two  models.  Figure  12  shows  the  resulting  electric  field  produced  at  a  9.75- 
mm  depth  for  both  models.  We  compared  the  electric  fields  at  31  other  representative 
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Electric  Field  (V/m) 


Time  (nsec) 

Fig.  12.  Electric  fields  for  two  models  of  water.  THe  fields,  which  result  from  the 
same  incident  pulse,  are  computed  at  a  9.75-mm  depth  for  each  model. 
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depths,  from  0  cm  through  4  cm,  and  observed  similar  agreement,  commensurate  with 
two  different  10%-accuracy  fits  to  the  water  data.  A  similar  degree  of  agreement  was 
also  seen  for  the  two  models  in  computations  for  the  following  incident  fields;  (1) 
A  1-cycle  8-GHz  square-modulated  sinusoid;  (2)  A  1-cycle  10-GHz  squcure-modulated 
sinusoid;  (3)  A  50-psec  square  pulse;  and  (4)  A  100-psec  squcire  pulse.  These  results 
axe  numerical  evidence  that  the  computed  interned  fields  are  stable  as  the  water  model 
is  perturbed  by  about  10%. 

Now  we  wish  to  provide  an  explanation,  based  on  the  understanding  devel- 
oped  herein,  of  the  observations  in  Ref.  3  whereas  no  significant  differences  in  SAR 
(Specific  Absorption  Rate)  distribution  between  pulsed  and  CW  (Continuous  Wave) 
expos\ires  were  measured  for  a  MEM  (Muscle  Equivalent  Material)  exhibiting  two  re¬ 
laxation  times.  The  dielectric  model  was  composed  of  two  Debye  mechanisms,  which 
fit  experimentally  determined  permittivity  and  conductivity  data  for  MEM  at  2.07, 
2.8,  5.6,  and  9.3  GHz.  The  two  relcixation  times  were  rj  =  6.63  psec,  and  T2  =  83.7 
psec.  The  material  was  illuminated  with  a  train  of  square-modiilated  pulses  of  var¬ 
ious  durations  and  repetition  rates.  We  will  explaiin  the  observations  in  Ref.  3  for 
the  smallest  pulse  repetition  rate  used  (200  pulses  per  second)  for  which  the  pulse 
duration  was  0.5  /rsec.  All  other  results  in  Ref.  3  with  different  pulse  settings  are 
similarly  explained.  For  this  incident  signal  the  author  of  Ref.  3  used  a  czirrier  of 
5.6  GHz.  Thus,  the  pulse  duration  was  Tp  =  5  x  lO”’^  sec,  the  quiescent  interval 
between  pulses  was  T,  =  5  x  10“®  sec,  and  each  pulse  contained  2800  cycles  of  car¬ 
rier.  Further,  it  happened  that  Tp  »  mas{r}  =  T2,  thus  the  medium  would  not 
respond  in  a  dispersive  manner,  rather  it  exhibited  an  effective  relative  permittivity 
of  ci  +  ej  —  foo  =  42.4  +  15.8  —  4.3  =  53.9  to  the  envelope  of  the  pulses.  With  these 
parameter  values  in  mind  we  expect  the  medium  to  sense  a  CW  signal  of  carrier  5.6 
GHz,  even  in  the  pulsed  case.  Since  the  MEM  has  low  heat  diffusivity,  i.e.,  it  takes 
about  40  sec  for  a  temperature  change  of  0.04  °C  to  occur,  its  temperature  will  change 
immezisurably  in  the  time  T,  between  pulses.  Consequently,  the  temperature  will  re¬ 
main  constant  until  the  next  pulse  arrives  again  to  be  seen  by  the  medium  as  peurt 
of  a  CW  exposure,  hence  the  observed  indistinguishability  of  the  CW  SAR  versus 
the  pulsed  SAR.  As  the  carrier  frequency  is  increased  the  CW  vs.  pulsed  exposure 
SARs  should  start  disagreeing  at  smaller  depths.  This  is  related  to  the  frequency- 
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dependent  skin  depth  for  the  carrier  component  which  is  9.7  mm  at  5.6  GHz.  All 
the  measurements  of  SAR  were  obtmned  at  depths  well  within  the  skin  depth.  The 
effects  of  the  pulsing  should  be  observable  at  depth  greater  than  C0T2  ~  1.2  cm  since 
then  the  Ccirrier  component  will  have  decayed  sufficiently  (it  is  exponenticiUy  decaying 
as  in  the  CW  case)  so  the  remaining  field  will  be  due  to  the  square  envelope  and  will 
behave  diffusively. 

The  results  presented  in  our  paper  also  help  in  accurately  predetermining  sam¬ 
ple  thickness  to  be  used  in  single  and  total  transmission  Time  Domaiin  Spectroscopy 
(TDS)  studies  such  cis  t^ose  presented  in  Ref.  26.  In  the  single  transmission  approach 
one  studies  the  first  arrival  through  a  long  sample  so  that  the  highest  frequency 
components  will  have  decayed  sufficiently  in  order  not  to  mask  the  lower  frequency 
components  which,  as  we  showed  in  Section  2,  are  significantly  slower  (for  water 
Cl  ~  Co/9).  The  first  suggestion  arising  from  the  aneilysis  of  Section  2  (verified  by 
the  numerical  simulations)  is  that  the  sample  can  be  as  short  as  2cor  when  one  is 
interested  in  measuring  the  static  permittivity  of  a  Debye-like  material  with  single 
transmission  TDS.  On  the  other  hand,  in  the  totad  transmission  approach  the  first 
zurival  time  through  a  short  sample  is  used  to  best  me2isure  the  infinite-frequency 
permittivity  of  the  material  under  test.  Therefore,  a  sample  length  shorter  than  CgT 
should  be  appropriate  in  order  to  capture  the  time  of  arrived  of  the  highest  frequencies 
and  thus  determine  Coo.  (The  next-to-last  peiragraph  of  Section  2  hcis  related  com¬ 
ments.)  For  the  test  case  presented  in  the  results  section  of  Ref.  26  [Eqn.  (8)  there] 
it  happens  that  e,  —  17.3,  £„  =  3.3,  r  =  460  psec.  Thus,  the  highest  frequency 
components  travel  with  a  speed  Cg  =  0.5505c,  the  lowest  frequency  components  travel 
with  speed  ci  =  0.2673c,  and  c^r  =  7.6  cm.  The  experimenters  could  have  used  a 
slab  of  material  of  at  least  15  cm  to  reliably  measure  and  a  slab  thickness  of  at 
most  7  cm  to  measure  £03,  instead  of  using  50  cm  and  2.7  cm  respectively  to  do  the 
measurements.  In  Figure  13  we  show  the  x-t  location  of  the  peak  of  the  impulse 
response  for  this  medium.  (The  two  data  points  in  one  corner  of  the  graph  are  minor 
numeriezd  artifacts  that  are  related  to  the  similcir  numerical  artifacts  in  Fig.  7.)  The 
graph  as  a  whole  confirms  our  predictions  of  the  sample  lengths  with  which  TDS  will 
be  most  successful  in  measuring  the  permittivities. 
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Kg.  13.  Th«  thne  of  airivol  of  the  of  the  Green  function  GE(r,t)  nt  various 
depths  in  the  medium  of  Ref.  26. 
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5.  CONCLUSION 

There  are  many  generalized  methods  for  computing  the  response  of  any  single  dis¬ 
persive  medium  to  any  single  incident  pulse.  This  paper  has  contributed  nothing 
along  those  lines;  instead,  we  have  developed  several  related  results  that  are  general 
in  a  different  way.  We  have  found,  for  instance,  that  the  power  density  and  the 
peak  amplitude  of  an  incident  pulse  place  upper  bounds  on  the  peak  amplitude  suf¬ 
fered  inside  a  dispersive  medium — independent  of  the  other  details  of  the  incident 

pulse.  Our  Green-function  results  reduce  the  computation  of  such  upper  bounds 

« 

to  little  more  than  a  calculator  exercise,  and  we  have  given  numerical  examples  in 
which  these  sharp  upper  bounds  are  almost  attained.  We  have  reported  similar 
upper-bound  estimates  for  a  quantity  related  to  power  density,  and  for  the  time 
derivative  of  the  magnetic  field.  That  time  derivative,  5tJT,  is  largely  responsible 

I 

for  electromotive-force  currents  in  circuit-like  structures;  it  is  especially  large  for 
short-risetime  pulses.  Such  upper  bounds  could,  potentially,  help  in  the  regulation 
of  electromagnetic  interference  or  damage  produced  in  dispersive  media. 

Although  our  methods  apply  to  dispersive  media  generedly,  we  have  used 
a  Debye  model  for  microwave-pulse  propagation  in  water  as  a  specific  numerical 
example.  For  that  water  example,  we  reported  a  factor-of-nine  effect  in  the  wave 
speeds  that  seems  to  have  been  unnoticed  until  now,  and  we  explained  this  large 
effect  analyticcdly  using  PDEs.  The  PDE  analysis  also  yielded  simple  short-pulse 
and  long-pulse  approximations,  as  did  an  ancdysis  of  the  numerical  Green’functions 
involved  in  the  upper-bound  concepts  described  earlier.  We  studied  rates  of  decay 
in  Subsection  2 A,  and  approximations  for  detailed  time  traces  are  given  in  Subsec¬ 
tion  2B  and  Eqns.  (3.2)  and  (3.4).  Section  4  uses  that  work  to  explain  the  results 
of  some  laboratory  experiments,  and  to  offer  suggestions  for  future  experiments. 
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APPENDIX  A:  ERROR  ESTIMATES 

This  appendix  develops  an  easily  used  method  for  estimating  the  percentage  error 
in  grid-dependent  computations.  The  estimates  are  validated  numerically  in  some 
cases  for  which  exact  solutions  are  known,  and  the  estimates  are  also  used  in  cases 
for  which  exact  solutions  are  not  known.  These  estimates  use  only  one  of  the  many 
definitions  of  relative  error  for  functions  of  two  variables,  but  the  estimates  can  easily 
be  adapted  to  other  measures  of  relative  error  and  to  functions  of  more  than  two 
variables.  We  were  motivated  to  estimate  percentage  errors  because  quantitative 
estimates  of  uncertainty  could  help  in  setting  safety  standards  or  in  the  use  of  other 
computations  for  which  measures  of  numerical  uncertainty  are  vital. 

The  error  estimates  in  this  paper  are  empirical.  The  estimates  are  obtadned 
by  computing  with  several  different  grid  sizes,  noticing  a  pattern  in  the  relative 
errors  of  the  different  computer  runs,  and  using  that  inferred  pattern  to  sum  the 
series  on  the  right  side  of  the  inequality  \\fc  -  /*  loo  <  1/c  -  /i  loo  +  l/i  -  /2  |oo  -r  I/2  - 
/3I00  +  •  •  •.  The  concepts  of  rate  and  order  of  convergence  are  not  used  in  these 
empirical  estimates.  The  logic  behind  the  estimates  is  slightly  intricate,  but  the 
resulting  method  uses  only  least  upper  bounds  and  geometric  means  and  geometric 
series,  which  are  extraordinarily  simple  concepts. 

We  will  estimate  relative  errors  for  functions  of  two  variables.  The  error  in 
a  computed  solution  relative  to  the  exact  solution  fc{z,t)  can  be  defined  as 

ll/c(^,‘)  -  /e(z,-)||oo 
i/c(^,-)|oo 

The  relative  error  (Al)  is  a  function  of  the  depth  z.  For  instance,  if  the  relative 
error  is  known  to  be  less  than  1%  over  a  range  of  depths  z,  then,  for  those  depths, 
the  exact  solution  fc{z,t)  can  differ  from  the  corriputed  solution  fc{z,t)  by  no  more 
than  1%  of  the  peak  value  (with  respect  to  f)  of  |/e(z,  t)l.  The  just-mentioned  peak 
value  of  |/e(z,-)l  is  unknown  in  all  practical  cases  in  which  the  exact  solution  /.  is 
unknown;  however,  in  the  present  hypothetical  case  of  1%  relative  error,  the  peak 
value  of  the  unknown  quantity  l/e(z,  •)]  can  differ  by  no  more  than  1%  from  the  peak 
value  of  the  known  quantity  l/c(z,  •)!.  A  brief  derivation  involving  a  geometric  series 
then  shows  that,  in  the  case  of  1%  relative  error,  the  graph  of  the  exact  solution 
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fe[z^t)  is  guaranteed  to  fall  between  the  graph  of  /c  (z,o  +  i.0I%||/42,-)I»  and 
the  graph  of  /c(z,t)  — 1.01%|/c(z,-)||oo-  The  quantities  ±1.01%|/c(z,  •)||oo  are  called 
error  bars. 

We  will  now  go  step  by  step  through  the  error  estimates  illustrated  in  Fig.  14. 
Those  estimates  are  for  Green-function  computations  of  waves  normally  incident  on 
water.  The  computations  were  done  for  five  different  grids,  and  the  results  are 
called  f2{z,t), . . . ,  in  order  of  the  increasing  fineness  of  the  grids. 

Because  of  computer  limitations,  the  finer-grid  computations  were  done  for  shal¬ 
lower  depths  than  were  the  coarser-grid  computations.  For  zero  depth,  the  light 
circles  on  the  2  =  0  axis  of  Fig.  14  represent  |/i(0,*)  —  /2(0, •)l|oo/|/5(0,*)|oo> 
1/2(0,-)  -  /3(0,-)1„/|/5(0,-)|„,  1/2(0,-)  -  /4(0.-)l„/||/5(0,-)l„,  and  l/,(0,-)  - 
/5(0,-)U/IIA(0,-)|oo,  running  from  top  to  bottom  along  the  vertical. axis.  For  in¬ 
stance,  the  error  in  relative  to  the  most  finely  computed  result  /5(0,t)  is 

0.157%.  The  even  spacing  of  the  four  light  circles  cdong  the  logarithmic,  vertical 
ajds  of  Fig.  14  suggests  that  the  four  relative  errors  are  in  geometric  progression; 
in  fact,  the  ratio  of  the  second-largest  relative  error  to  the  largest  relative  error 
is  0.368,  the  ratio  of  the  next  two  smzdler  relative  errors  is  0.284,  the  ratio  of  the 
two  smallest  relative  errors  is  0.268,  and  the  geometric  mean  of  the  three  previous 
numbers  is  0.304.  For  those  reasons  we  assume 


l/n-n(0>‘)  /n-4-2(0}  Olloo  ^  ^041"  l/l(Q>*)  /2(0|*)l|oo 

11/5(0,.)  loo  1/5(0,  •)  Hoc  ’  ^  ^ 

even  for  hypothetical  results,  such  as  /loo,  which  would  come  from  computations 
involving  a  much-finer  grid  than  was  used  for  the  actual  computation  of  We 
also  assume  that  the  exact  result  /e(2,t)  is  the  n— >oo  limit  of  /„(z,t).  The  triangle 
inequality  for  norms  then  implies  |/i  -/eioo  <  l/i  -/2IC0  +  I/2  -/3ioo  +  |/3-/4ioo  + 
• . ..  We  use  (A2)  in  the  right  side  of  the  previous  inequality  to  get  a  geometric  series, 
which  sums  to 


ll/i(0,.)-/c(0,' 


II 00  < 


1 


11/5(0,  .)!< 


1  -  0.304 


II/i(0,-)-/2(Q.' 
1/5(0,  •)||co 


=  8.08% 


(.43) 
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at  z  =  0.  That  8.08%  estimated  relative  error  is  plotted  as  the  top-most  heavy  circle 
on  the  verticzJ  axis  of  Fig.  14.  The  analysis  was  repeated  for  the  relative  errors  of 
other  com.putations  fi,  yielding 

Il/.(0,.)U - - 0i(MIU - ’ 

for  i  =  1,2, 3,4.  Those  estimated  relative  errors  are  also  plotted  as  heavy  circles  on 
the  vertical  axis  of  Fig.  ^14. 

Reference  9  has  an  exact,  closed  form,  modified-Bessel-function  expression 
of  the  zero-depth  Green  function  for  a  hadf-space  of  Debye  medium.  That  exact 
solution  /c  Wcis  used  to  compute  the  true  relative  errors  Erc\[fi,fe]  at  z  =  0.  Those 
true  relative  errors  are  plotted  as  the  bold  Xs  on  Fig.  14.  The  true  relative  errors 
(bold  Xs)  validate  the  error  estimates  (bold  circles)  because  the  two  sets  of  errors 
match  closely;  for  instance,  the  error  estimate  (A4)  for  /3(0,  •)  was  that  there  would 
be  0.844%  error  relative  to  /*,  and  the  actual  error  in  /3(0,-)  relative  to  /e  was 
0.797%.  Our  percentage-error  method  involving  geometric  means  and  geometric 
series  is  thereby  validated. 

The  relative  errors  lfi{z  r)-  fi+l{^,-)loo/lfs{z,-)\\oo  were  then  computed 
at  the  depths  z  =  0.480  mm,  1.04  mm,  1.52  mm,  and  2.00  mm,  for  i  =  1,2, 3, 4. 
Those  results  are  plotted  as  some  of  the  light  circles  to  the  right  of  the  vertical 
axis  in  Fig.  14.  Error  estimates  for  fz{z,‘)  were  computed  using  the  geometric- 
series  technique  in  (A2)-(A4),  with  the  results  plotted  as  severed  bold  asterisks  on 
Fig.  14.  The  ratios  of  relative  errors  and  the  geometric  means  of  the  ratios  were 
computed  separately  for  each  depth.  In  particular,  the  geometric  means  for  the  four 
depths  mentioned  in  this  paragraph  are  0.302,  0.328,  0.315,  and  0.333  for  0.480  mm 
through  2.00  mm,  consecutively.  We  will  now  explain  how  the  relative  errors  for 
depths  beyond  2.00  mm  were  computed. 

The  finest-discretization  run  was  not  computed  beyond  2.00  mm  for  reasons 
related  to  computer  resources.  Thus,  only  /i,  fz,  fz-,  and  were  anadyzed  for  the 
depths  from  2.48  mm  through  4.00  mm.  These  analyses  were  done  as  described  in 
the  previous  paragraph,  and  the  geometric  means  were  adso  computed  separately 
for  each  depth.  Similarly,  the  runs  /i,  fzt  and  fz  were  anadyzed  for  depths  from 
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Relative  Error  (1  =  100%) 


Depth  (millimeter) 

Fig.  14.  Relative  errors  for  the  computation  of  Gb  for  water.  The  bold  marks  rep¬ 
resent  the  final  error  estimates  and  their  numerical  validations.  The  lighter 
marks  represent  intermediate  steps  in  the  computation  of  error  estimates. 


0 


5.04  mm  through  7.04  mm.  In  these  computations,  and  in  all  other  computations 
described  in  this  appendix,  the  Loo  norms  were  computed  using  90  evenly  spaced 
time  points. 

The  final  result  in  Fig.  14  is  symbolized  by  the  bold  asterisks  there.  The 
result  is  that  fz{z,t)  has  errors  relative  to  that  are  expected  to  fall  monotonically 
from  about  1.70%  at  z  =0.480  mm  to  about  0.799%  at  z  =4.00  mm.  The  error  in 
fz{0,t)  relative  to  /e(0,t)  is  actually  0.797%.  Figure  14  shows  that  the  relative 
errors  are  non-monotonij  from  0.00  mm  through  0.500  mm;  therefore,  we  make  no 
further  inferences  about  the  relative  errors  in  that  first  half-miUimeter. 
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